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DEFINITION OF SYMBOLS 


English Letters 


A 


B 


Area 


A ,A _ 
P’ PB 




area of convection; subscript B denotes value 

at the updraft base level 

area of subsidence or cloud-free air 

area represented by a grid point, l.e. Ax»Ay 


Bowen ratio 


Buoy parcel buoyancy as defined on p. 41 

C linear advectlve wave phase speed 


convective parcel aerodynamic drag coefficient as dis- 
cussed on p. 44 


CG 


specific cloud water growth rate as defined on p. 8 


c 


c 

8 


c 

V 


specific cloud or suspended water (mass per unit mass of 
moist air) 

c , Kessler conversion threshold value 
th 

gravity wave phase speed 

specific’ heat of dry air at constant pressure and con- 
stant volume, respectively 


pm 


specific heat of moist air at constant pressure 


£ 


transition layer height 




vapor pressure; subscript s denotes saturated value 


f 


Coriolis parameter 


f 

CO 


fraction of updraft core without rainwater falling 
through it 




Integrated form of the Buslnger relationships for heat, 
momentum and specific humidity, respectively 


f 


ro 

g 

H 

«s»»so 


fraction of rainwater falling outside the updraft core 
acceleration due to gravity 

mean height used In linear shallow-water equations 

surface available energy flux as discussed on p. 13, sub- 
script 0 denotes value at the equator 
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h 


surface layer depth 


h 


S 


K 


K 

z 




k 


sensible heat flux through level z 


0 


conv 

coll 

evap 

fall 


Kessler cloud physics parameterization coefficients 
conversion coefficient used in Chapter 2 
collection coefficient used in Chapter 2 
evaporation coefficient used in Chapter 2 
modified terminal velocity coefficient used 

in Chapter 2 

Kj, conversion coefficient used in Chapter 3 

K 2 modified terminal velocity coefficient used in 
the collection calculation in Chapter 3 
collection coefficient used in Chapter 3 
K, modified rainwater- fallout coefficient used 

in Chapter 3 

modified rainwater evaporation coefficient 
used in Chapter 3 

K, rainwater evaporation coefficient used in 

Chapter 3 


vertical 

^z 



component of the diffusion coefficient 
vertical component of the diffusion coefficient 
for heat 

vertical component of the diffusion coefficient 
for momentum 

vertical component of the diffusion coefficient 
for specific humidity 


horizontal components of the filter coefficient for the 
fourth- and second-derivative filters 


von Karman constant (0.35) 


L Monln-Obukhov length 

L ,L^,L^,L^ latent heat of condensation, fusion and sublimation, 

respectively. L is either L or L depending on the 
temperature 


vol 


characteristic length of the plume parcel (cube root of 
the parcel volume) 


M 


P 


mass 




^OP 


"fuel" mass which feeds the plume 
plume mass which passed through a level; sub- 
script B denotes value at updraft base level 
plume mass which passed through z,pQp 


total plume mass in a grid volume 


X 



total subsiding mass In a grid volvune 
specific humidity flux through Zq 
number of clouds In a grid volume 

number concentration of rain drops as discussed on p. 62 

precipitation amount at the updraft base level and at the 
surface, respectively 

convective parcel depth 

pressure 

p^ pressure at the model top, 

p^ subsided environmental pressure 

nonconvectlve dlabatlc heating rate 

specific humidity growth rate as defined on p. 8 

specific humidity Cinass per unit mass of moist air) 
q specific humidity of the plume environment 

specific humidity at level h 
q specific humidity of the plume 

q^ saturated specific humidity 

qg surface specific humidity 

q saturated specific humidity of the plume environ- 

S6 . 

raent 

q^p saturated specific htjmidity of the plume 

Convective cloud or suspended water concentration (mass 

per unit mass of moist air) 

q cloud water in the plume environment 

q cloud water in the plume 

^cp 

q Kessler cloud water conversion threshold value 

*^th 

convective precipitation water concentration dropped or 
unloaded by the ascending parcel (mass per unit mass of 
moist air) 

convective rain or precipitation water concentration in 
the ascending parcel (mass per unit mass of moist air) 

convective ice concentration in the plume (mass per unit 
mass of moist air) 

convective rainwater concentration in the plume "aging" 
process O^ass per unit mass of moist air) 


<1* 

R 



R 


V 


RG 

RH 

r 



t 


U 

u 


V 


W(I),W(x) 


surface friction specific humidity 
radius 

R^ updraft base radius 

R^ initial updraft base radius 

R^ radius of influence for cloud 
R^qP cloud radius at the cloud top, z^Qp 

specific gas constant for dry air and water vapor, 
respectively 

rainwater growth rate as defined on p. 8 
relative humidity 

specific rain or precipitation water (mass per unit mass 
of moist air) 

radius of the earth 


layer releasable instability as defined on p. 42 
temperature 

Tjj final freezing temperature 

initial -freezing temperature 
T moist or wet temperature as defined on p. 7 
T subsided environmental temperature 
T° virtual temperature 

environmental virtual temperature 
T parcel virtual temperature 

subsided environmental virtual temperature 

time 

local time 

tj time constant for the parcel freezing rate 
mean velocity used in linear shallow-water equations 
west-to-east component of the horizontal wind 
surface friction velocity 

horizontal vector wind 

horizontal vector wind at level h 
n 

Vg surface horizontal vector wind 

terminal fall velocity of precipitation 

lateral boundary discrete and analytic weighting 
coefficients 


xii 



I 


vertical velocity 


,w ,w 
cs’ cu 


w ,w 


P pi 
pmin 


mean cumulus-scale vertical velocity, mean 
cumulus-scale subsidence rate and mean 
cumulus-scale upward motion as discussed 
on p. 52 

convective parcel vertical velocity; sub- 
script 1 denotes Initial Impulse velocity 
minimum parcel vertical velocity 
net convective scale vertical velocity as 
discussed on p. 64 


X 

y 

z 


the west-east horizontal coordinate 

the south-north horizontal coordinate 

height above the ground and the vertical coordinate 


“B 

'top 


height of the cloud base 
height of the cloud top 
height of the FE model top 


surface roughness length 


Greek Letters 


a 


a 


P pB 


a 


r 

Y 

^e 

A 


convective parcel entrainment parameter 
earth's latitude 

percent area of cloud; subscript B denotes value at 
updraft base level 

percent area of subsiding or cloud-free air 

filter function stability coefficient 
(4) 

B filter function stability coefficient for.fourth- 
derivative filter 

B'' filter function stability coefficient for second- 
derivative filter 

adiabatic virtual temperature lapse rate 
virtual mass parameter as discussed on p. 44 
environmental virtual temperature lapse rate 
grid interval 

Ax west-east grid Interval 
Ay south-north grid Interval 
Az vertical grid Interval 


xiii 


Increment of cohvectlve rainwater in the plume "aging" 
process used for collection and evaporation calculations 

difference between the saturation specific humidity over 
water and over ice 

time step 

vertical hydrostatic pressure difference in the plume 
between Zg and 

vertical hydros^tic pressure difference in the subsided 
environment between z and z„_„ 

D roiT 

depth of the convective fuel layer 

subsidence in cloud-free air 

number of Aq^ bundles in q 

potential temperature 

6^ potential temperature at level h 
_0g surface potential temperature 
0g mean surface layer potential temperature 
0^ virtual potential temperature 

surface friction temperature 

wavenumber 


convective parcel entrainment rate 

time-smoothing coefficient as discussed on p. 23 

calling frequency of the convection parameterization as 
discussed on p. 67 

ratio of the radius of a circle to its circumference 


virtual density as defined on p. 7 
density of dry air 
virtual density of the environment 
moist or wet density as defined on p. 7 
moist density of the plume at the updraft base 
the plume environment 
the plume 
the plume at z 
the subsided environment 
of the subsided environment 


3 ™ 

mB 

) 

^me 

mp 

^mpTOP 

ms 


moist density of 
moist density of 
moist density of 
moist density of 
virtual density 


w 


density of liquid water 


cumulus time scale 



T 


s 


momentiun flux through 


z 


0 




normalized gradients of heat, momentum and specific 
humidity 


X 


symbol used to represent a group of variables as defined 
when used 


b) 


dp/dt 


Subscripts 


c 




( 



variable evaluated at grid point l,j,k or position 
Xj^+(l-l)Ax, y^+(j-l)Ay, Zj^+(k-l)Az where x^, y^^, Is 
the origin 

variable evaluated at vertical grid point corresponding 
to the model top, z^ 


Superscripts 

'I 

( ) variable evaluated at time step T 
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1.0 INTRODUCTION 


1.1 Statement of the Problem 

The prediction of convective activity and Its severity has been a 
problem for many years. These forecasts are needed so the public can 
prepare for the possible disastrous consequences of these storms. The 
draft of the NOAA Project Development Plan for Project Revere E^nviron- 
mental ^torms jmd Mesoscale Experiment (SESAME) (1974) states; "In the 
7-year period from July 1963 to June 1970, for example, severe local 
storms In the United States alone accounted for over 2300 deaths and 
4 billion dollars damage." 

It is possible to divide convective activity into two categories 
based on the physical mechanism initiating the activity; 1) thermally 
forced convection and 2) dynamically forced convection. This categori- 
zation is not meant to imply that these forcing mechanisms are mutually 
exclusive; it is clear that in nature there is a continuous transition 
from one category to the other. 

Examples of thermally initiated convection include afternoon air- 
mass convection, convection caused by cold air moving over warm water, 
i.e. cases where the surface heat flux destabiizes the environment such 
that surface rooted convection can occur. In particular, many of the 
hail- producing convective cells over the National Hail Research 
meat (NHRE) region in Northeast Colorado occur late in the afternoon 
after the surface heat flux has increased the depth of the well-mixed 
layer to the extent that it coincides or at least nearly coincides with 
the updraft parcel’s lifted condensation level, LCL (Perkey, 1973). 
Until the surface adiabatic layer builds high enough to encompass the 
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LCL, little, if any, significant convection can occur. The thermal- 
induced activity can range in severity from innocent ^afternoon showers 
to giant hail -producing thunderstorms which can cause crop and property 
damage. 

Examples of the dynamically forced convection include the frontally 
forced squall line (Palmen and Newton, 1969), the sea-breeze convergence 
convective band (Pielke, 1973) and warm frontal convection embedded in 
altostratus (Kreitzberg and Brown, 1970). In these cases, the lifting 
due to the frontal convergence converts the potential instability to 
realizable conditional instability which results in convection. The 
severity of these storms runs from the mild afternoon sea-breeze shower 
to squall-line tornadoes and flood-producing warm-frontal convection. 

Observations show that mesoscale organization of dynamically forced 

convective storms often occurs within cyclonic scale systems. These 

mesoscale areas or bands of showers have horizontal dimensions typically 

50 to 200 km by 500 to 800 km and are characterized by convergences on 
-4 -1 

the order of 10 s , deeper moist layers than non-convectlve regions, 
extensive cloudiness and embedded clusters of showers (Elliott and 
HoVind, 1964; Kreitzberg, 1964; Kreitzberg and Brown, 1970; Browning arid 
Harrold, 1969; Browning ^ al. , 1973; Matsumoto, Nlnomiya and Akiyama, 
1967). 

Observations indicate, therefore, that the convection oftentimes 
builds in a mesoscale environment with different thermal properties and 
more low-level convergence and hydrostatic vertical motion than is 
characteristic of the cyclonic scale. One, therefore, is faced with 
cyclonic scale systems with 1000 km wavelengths and vertical motions of 
a few cm/ s in which are embedded essentially hydrostatic mesoscale 
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systems with wavelengths of 100 km and vertical motions of a few tens 
of cm/s. Within the mesoscale disturbed areas are embedded non- 
hydrostatlc convective systems with wavelengths of 10 km and vertical 
motions of a few m/s. 

At this time, it is not practical or possible to explicitly calcu- 
late details of the convective activity; it must be parameterized. In 
order to predict the range of convective activity which includes both 
the thermally and dynamically forced convection, this parameterization 
scheme must be general. Also, the driving model must have a fine 
enough grid to resolve length scales on the order of 200 km up through 
the synoptic scale, thus yielding convergence and vertical motion 
patterns with enough detail to represent the environment in which the 
convection must evolve. 

1.2 Research Obiectlves 

This report concentrates on the dynamically forced convective 
activity and, in particular, on the frontally forced squall-line pre- 
cipitation region. It is the squall line which displays the most 
violent behavior and, thus, the most potential for damage. These are 
the storms which can produce locally heavy rains and their associated 
flash floods, which can bear hall that produces crop and roof damage, 
and which can spawn tornadoes and their destructive winds. 

To attack this problem, a parasitic-nested, fine-mesh, three- 
dimensional, moist, primitive equation model has been developed. This 
model contains a general convective parameterization scheme which has 
been developed such as to be applicable for both the thermally and 
dynamically forced convection. The details of this scheme and how It 
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differs from previous schemes Is discussed In Chapter 3. Because of Che 
local nature of the squall line, a mesoscale grid (A " 35 km) is desir- 
able; however, for the tests in this report, a fine-mesh grid 
(A •” 132 km) is adequate. Even the fine mesh necessitates the use of 
a limited-area model which must reflect the time changes on the large 
scale. To satisfy this need, a set of time-dependent lateral boundary 
conditions has been developed and tested. See Section 2.3.3 and 
Appendix 2 for a discussion of these conditions and their relation to 
other boundary conditions presently in use. 

Using the real-data case of 21 February 1971 this report 
addresses several questions: 

1. Are the lateral boundary conditions as proposed acceptable for 
real-data forecast use? 

2. Is the convective parameterization scheme general for mid-latitude 
short-term forecasting; l.e., will it yield reasonable results in 
both the frontally forced squall-line region and the higher based 
stratus embedded convective regions near the center of the low? 

3. What is the effect of the horizontal observational resolution on 
the 12 h to 18 h forecast quantitative precipitation amounts? 

4. What is the effect of the model vertical resolution on the 12 h to 
18 h forecast quantitative precipitation amounts? 
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2.0 THE PRIMITIVE EQUATION MODEL 

Several choices must be made before and during the development of 
a numerical model. Many of these choices must be made arbitrarily, while 
others can be based upon the experience and results of other researchers. 
Thus, while all numerical models have much in common, there are numerous 
differences which affect their behavior and suitability for specific 
meteorological problems. 

The basic characteristics of this model are: 

1. The height above a spherical surface is used as the vertical co- 
ordinate. Terrain effects are not included. 

2. The model is hydrostatic; thus, the horizontal winds (u,v), virtual 

temperature (T^) , specific humidity (q), cloud or suspended water 

(c), and precipitation water (r), are predicted. In addition, 

pressure at the model top (p ) is predicted assuming u) =0. 

^t ^t 

Pressure at all other levels is diagnosed hydrostatically. Density 
(p) is diagnosed using the gas law. Vertical velocity (w) is diag- 
nosed from Richardson's form of the continuity equation. The 
effects of water loading are included in the hydrostatic calcula- 
tions. 

3. The convective transports are evaluated using an adjustment scheme 
which is based on a one-dimensional Lagranglan cumulus model. 

4. As mentioned above, the upper boundary condition is (dp/dt) = 

^t 

0 ) = 0 . 

^t 

5. The lateral boundaries use a blending of specified and model cal- 
culated tendencies on all prognostic variables. Also, a boundary 
region of high viscosity is needed. 
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6‘. The lower boundary employs a Monin-Obukhov (1954) surface layer 
with an O'Brien (1970) K^-profile in the transition layer. 

7. The horizontal filtering necessary in numerical models is calcu- 
lated using a fourth-derivative diffusive filter. 

8. The finite-difference formulation is basically fourth-order 
centered in space and leapfrog in time. A weak time filter is 
used to avoid separation of solutions. 

2.1 The Primitive Equations 

The primitive equations are cast in an x,y,z coordinate system on 
a latitude-longitude grid with first-order corrections for the spheri- 
cal shape of the earth. The effects of terrain are neglected throughout 
the calculations. The momentum equations are 


3t 


3v 

3t 


- w + fv 


1 9p . uv ^ \ , 

— — ^ tan (ot ) + 

p 3x r e 

m e 


3u \ 

3t/ ' 

' c 


^ 3v ^ 1 

-V»Vv - w X fu 

3z p 


m 


3p U fr, \ J. 

^ tan (a ) + v- 

3y r e ot 




( 1 ) 

( 2 ) 


where r denotes the earth's radius, a denotes the latitude and the s 
e e 

subscripted tendencies refer to sub-grid scale tendencies. The p^ is 
the wet or moist density and includes the effect of liquid water load- 
ing (see definition below). The vector notation refers to the hori- 
zontal components only. The thermodynamic equation is written as 


3T 
V 

3t 





3T, 
\ 

3t 


conv 



( 3 ) 
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where 


P 


R.T ’ 

d V 


= T(l+0.61q), 


m 


R.T 
d m 


T 

m 


T 

V 


1+c+r 


. » 


m 

Q is the nonconvective diabetic heating rate, c the cloud water and 
r the precipitating water. The local pressure tendency used in Eq. (3) 
is 


1e 

3t 


gP^w - g 




P V + - — tan(a ) 
m r e 


P ^ (c+r) 


dz 


*^p), 


(4) 


where oi =0.0 and z is the top of the model. For a derivation of 
temperature and pressure tendency equations with liquid water loading 
terms, see Appendix 1. At z^, the pressure is forecast using 


3 


l) ■ + SP„“), 




conv , z 


-I?) 


s.z. 


(5) 


The moisture equations are 


|a = _^.^q _ w |a + QG + + If ) . (6) 

'conv 's 

|£ . - „ || + CG + If ) , (7) 

'conv 

|| . 4.?, .„||+RG . (8) 


The conv subscripted tendencies on the right-hand side of the virtual 
temperature, pressure and moisture equations denote the changes in the 
quantity due to convection. The moisture growth terms QG, CG and RG 



8 


are given by 


QC 

CG 

RG 


- condensation rate + evaporation rate of c + evaporation 
rate of r, 

condensation rate - evaporation rate of c - conversion 
rate - collection rate. 


conversion rate + collection rate - evaporation rate 


of r - 


1 _ ^ 

p 8z 
m 


m t 


» 


where is the terminal fall velocity of the precipitation and is 
positive downward. 


V 


t 


^fall 


0.125 

r 


> 


where 

5.094 , T ^ 273. 15K; 

5.094[0.2+0.12*(r-0.5)] , T< 273. 15K . 


fall 


is the value suggested by Kessler (1969) when the rain is liquid 
but has been modified to account for the slow fall velocities of snow. 

A more complete discussion of this is in Section 3.5. 

The collection, conversion and evaporation rates are parameterized 
after Kessler: 


Conversion rate 
Collection rate 


K (c - c , ) , 

conv th 


^ .875 

coll ’ 


K (q - q)r 
evap ^s ^ 


.65 


> 


Evaporation rate of r 
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where and are prescribed following Kessler. 

The condensation rate and evaporation rate of c are calculated to main- 


tain saturation when suspended water is present and to avoid super- 
saturation. 

The diagnostic continuity equation is 


^ _2 V_ 

3z ^ cT cp 


fvi 

"pP 


P V 
z 


(V*^p) - (V»^p), 


/ t Pn.' 

5-p^v - — 


tan(a^) ~ P ^ (c+r) 


dz + — tan (a ) 
r e 

e 


(9) 


and the hydrostatic equation is 


^ = - -2— , (10) 

R^T 
d m 

where T^ is the layer-average T^. 

2.2 The Planetary Boundary Layer Formulation 

The planetary boundary layer is defined as the layer in which sur- 
face frictional forces cause the flow to deviate from gradient flow. 

For the purposes of modeling, this layer is divided into a shallow sur- 
face layer of depth h in which the fluxes are assumed constant and a 
transition layer of depth, E-h, in which the fluxes decrease in magni- 
tude to their f ree-atmosphere values (see Fig. 1). The depth of the 
constant flux layer, h, is somewhat arbitrary; however, if the selected 
h is too large, problems can occur at night when the atmosphere becomes 
quite stable. This problem has been discussed by Webb (1970), Taylor 
(1971) and others. Physically the trouble arises when L, the 



25 m 


1.22 m 


level 2 (z 2 “h) 


level 1 (z^) 


0.3 m 


0.0 


level S (Zq) 


7"7"7T7TT7"7TTT7TTT7 


GROUND 


Figure 1. Vertical grid for surface layer as used for the 
forecasts discussed in Chapter 4. 


Monin-Obukhov length, becomes much less than h. This implies that the 
boundary layer air and the air at level h are not coupled through mixing. 
Thus, the Buslnger log-linear relationships (Businger e^ al. , 1971) are 
no longer appropriate for diagnosing the meteorological quantities at 
1.22 m from the quantities at h. Webb (1970) suggested using a linear 
relationship when z/L is greater than one; however, the same difficulty 
arises when z/L becomes much greater than one. 

For the forecasts discussed in Chapter 4, h is set to 25 m. This 
is still too large for the stable night time regime but it does not 
cause severe enough problems to negate the forecast’s value. After 
Inspecting several 1200 GMT February 21 and 0000 GMT February 22 sound- 
ings for the southeastern U.S., the mixed layer height, E, was choosen 
as 1.0 km. A possible Improvement would be to Include a prognostic equa- 
tion for E £ Deardorff (1972). This approach would allow a diurnal 

; 1 f 

variation of E. The roughness length was set to 30 cm (Hoxit, 1974). 
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The general boundary layer formulation is similar to that of 
Pinkerton (Pinkerton and Kreitzberg, 1973; Kreltzberg et al. , 1974). It 
is assumed that K-theory provides an adequate and appropriate vehicle 
for calculating the vertical sub-grid tendencies included in the PE 
model, i.e. it is assumed that the vertical divergence of the time- 


averaged covariances such as 9u*^/9z, 9'^^'^/9z, 9q0^9z, can be re- 
presented by |j), II), 1^ II). The following 

notation has been used above and will be adhered to throughout the 
remainder of the discussion: the prime refers to scales of motion too 

small to be explicitly resolved by the model; the hat refers to the 
value of the quantity averaged over one grid-volume and one time step; 
the subscript M denotes momentum, the H denotes heat, the Q denotes 
moisture and the z denotes the vertical component of K. 

2.2.1 The Surface Layer 

The surface layer, z^ (the roughness length) to h, is modeled 
following Monin-Obukhov (1954) similarity theory, i.e, the surface 
gradients are written as follows: 



(lla) 

(llb) 

(llc) 


where k is von Karman’s constant and u^^, 0^ and q^ are the scaling 
velocity, temperature and specific humidity, respectively. These 
scaling quantities are related to the surface fluxes of momentum, 
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sensible heat and moisture as noted: 


P“*» hg = -c u^0*P: 


S “S p 

lUg = -u^q^. (J)^, <|)jj and (p^ are non-dimensional functions related to 

the atmospheric stability; thus, accounting for buoyancy effects, L is 
the Monin-Obukhov length which is defined as 


kghe 


( 12 ) 


where hg is the heat flux through level and 6 g is the mean surface 
layer potential temperature. At a height z^, the surface horizontal 
velocity, V_, is assumed zero while the surface potential temperature, 
9g, and specific humidity, qg, are diagnosed. 

Buslnger et al. (1971) determined that the (}>'s should have the 
following forms, for stable conditions: 


*«(t) - ’ 

"'0-74 + 4.7(f) : 


fpr unstable conditions: 



(13a) 

(13b) 

(14a) 

(14b) 


(J>^ is assumed equal to <j>„. The surface stability is determined by the 

g H 

sign of the surface heat flux; positive flux is unstable, negative is 


stable. 
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X3 

As a first cut at calculating hg, the total surface available 
energy flux, H„, Is Introduced. This quantity Is the energy flux avail- 

D 

able from the surface to modify the air. In a complex surface radiation 
and energy budget calculation, Hg Is the resultant of the balance of the 
solar energy flux, the long-wave energy flux and the soil energy absorp- 
tion flux. For the calculations which follow, the total surface avail- 
able energy flux Is calculated analytically as the following function of 
local time and latitude; 

where t^ Is the local time and Hg^ Is the available energy flux at the 
equator. contains the ft?llowlng latitudinal variation; 

for 0000 ^ 1 0600 , 

for 0600 < < 1800 , 

for 1800 ^ < 2400 , 

where Is the grid point latitude. For the February case discussed 
In Chapter 4, Hg^ Is as follows; 

{ 1 j (m^ s)~^ for 0000 ^ t^ ^ 0600 , 

15 j (m^ s)”^ for 0600 < t^ < 1800 , 

1 j (m^ s)”^ for 1800 ^ < 2400 . 



The above daytime value Is approximately 10% of the average energy flux 
during February as calculated by mean energy budget methods. However, 



because of the extensive overcast cloud cover over the central and 
eastern U.S., the reduction of the available surface energy seems 
reasonable. 

When H Is positive, the available flux is split into sensible and 
latent flux using the Bowen ratio which is set constant in both time 
and space; B = 0.74. This value, according to Sellers (1965, p. 105), 
is the North American average. The sensible heat flux is calculated as 
follows 

" “s ( 1^ ) ‘ 

When Hg is negative, hg = Hg. 

The specification of the surface energy flux independent of the 
lower layer stability can cause the diagnosed surface temperature to be 
unrealistic in order to support the imposed flux. This Inconsistency 
occurs most frequently at night under very stable regimes. It seems 
physically possible to impose any surface flux desired in a closed-box 
experiment; however, the resulting profiles may not be log-linear as 
Buslnger’s empirical relationships imply is the case in the atmosphere. 
It should also be noted that often the Imposed negative heat flux causes 
a breakdown in the diagnostic relationships due to z/L becoming too 
large as discussed earlier. To overcome these difficulties, a lower 
limit is placed on L; L cannot be less than h. 

Eqs. (11) are Integrated from to h to yield 






Wing section drag coefficient, 



Wing-Body Terms 


, = 0.00639 

Cp, = -0.00013 

wing 

wing-on-body 


interference 

, = 0.00072 

Cp, =0.00013 

^body-on-wing 

body 


interference 


Nacelle Terms 



Inboard 

Outboard 

Isolated 

0.00075 

0.00075 

wave 

Body-on-nacelle interference 

-0.00002 

0.00000 

Nacelle-on-body interference 

0.00005 

0.00010 

Nacelle-on-nacelle interference 

Direct 

0.00034 

0.00023 

Image 

0.00054 

0.00046 

Wing-on-nacelle interference 

-0.00043 

-0.00058 

Nacelle-on-wing interference 

-0.00156 


I Cp) = 0.00064 
nac 


2 Wing-body-nacelle Crj =0,00775 
^wave 


FIGURE 3.6-2.-TYPICAL WAVE DRAG COEFFICIENT SUMMARY 
NEAR-FIELD PROGRAM (M = 1. 1) 


analysis programs. (This option is described in section 
3.7, but basically requires that the total surface 
pressure coefficient on the wing, i.e. , thickness+lift, 
cannot be less than some specified fraction of vacuum 
pressure coefficient.) 

If the vjing thickness pressures are to be used by the 
wing design or lift analysis programs in pressure 
liiriting options, then the near-field program must first 
be run. Curing program execution, the thickness 
pressures are loaded into a system common block and are 
then available where needed. 

Nac§ll§ p ressur e field option s . - The near- field program allows 
for up to 3 pairs of nacelles located external to the wing- 
fuselage (or 2 pairs plus a single nacelle at Y=0) . The nacelles 
may be either above or below the wing (or both) . 

The nacelle pressure field is the pressure field imposed on the 
surface of the wing by the nacelles. A feature of the near-field 
program is the choice of *'v7rap” or "glance" solutions for the 
nacelle pressure field, as shown in figure 3.6-3. (The far-field 
wave drag program uses essenrially the ’’wrap" solution) . 

Available experimental data do not make it clear whether a "v/rap" 
or "glance" solution is more correct. Since the nacelle-on-vring 
interference term is substantial, both solutions are available in 
the program (controlled by an input code) . 


3.7 Wing Design and Lift Analysis 

The wing design and lift analysis programs are separate lifting 
surface methods which solve the direct or inverse problem, of: 

« Design - to define the wing camber surface shape 

required to produce a selected lifting pressure 
distribution. The wing design program includes methods 
for defining an optim.um pressure distribution. 

9 Lift analysis - to define the lifting pressure 
distribution acting on a given wing camber surface 
shape, and calculate the associated force coefficients. 

The lift analysis program contains solutions for the effect of 
fuselage, nacelles, canard and/or horizontal tail, and wing 
trailing edge flaps or incremental wing twist. Using 

superposition, the program solves for drag-due-to-lift , lift curve 
slope, and pitching m.omenr cnaract eristics of a given 
conf iguration through a range of angles of attack at a selected 
Mach number. 
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PRESSURES "GLANCE" AWAY FROM WING AT ADJACENT 
NACELLES 



The nacelle pressure field and accompanyinp shock waves "glance" 
away from the wing when encountering adjacent nacelles. In 
application, the nacelle generated pressure field is terminated on 
encountering another nacelle. 



PRESSURES "WRAP" AROUND ADJACENT NACELLE 

The nacelle pressure fields and accompanying shock waves 
"wrap" around adjacent nacelles. In application, the nacelle 
generated pressure field is allowed to pass through another 
nacelle as if it were transparent. 


FIGURE 3.6-3 -NACELLE PRESSURE FIELD CONCEPTS 


19 

Energy amplification (\^) ■ Present energy (X ) ^ 

Initial energy 

2 # 

where the energy is defined as x )/2. This quantity is calculated 
for each wavenumber, X . The response curves for both the second- and 
fourth-derivative filters as used are shown in Fig. 2. It should be 
noted that the 3^^^ for grid point B+6 is chosen as the average of the 
values at grid points B+5 and B+7. 

2.3,2 The Finite-Difference Formulation 

2. 3. 2.1 Differencing for all variables except rainwater 
The horizontal and vertical advectlve terms for all variables except 
r are approximated by centered finite-difference schemes. The horizontal 
derivatives employ a fourth-order accurate representation, for example. 


lx 

dx 






At the boundaries, a one-sided fourth-order scheme is applied; the 
following is an example of the scheme as applied to the west boundary, 
B: 


lx 

3x 




, j "^^^B+2 , j 3 , J 4 , J 


]/12Ax^ , 


and applied to the grid points one Ax east of the west boundary, B+1: 


lx 

3x 


w-i.j 


j“10XB+i, j+lBXB+2,j“^XB+3 j+Xb+, 4, j ]/12Axj 


The vertical derivatives are approximated by a centered second- 


order accurate method for an uneven grid; 




to 

o 




(a) GRID NUMBER (b) WAVENUMBER 


Figure 2. Response curves for fourth-derivative filter with 3^^) = 0.005 (used for the interior 
grid points); fourth-derivative filter with = 0.060 (used for grid points B+2 

through B+5); and second-derivative filter with 3^^^ = 0.24 (used at grid .point B+1); 
(a) X (T'j) energy amplification. 
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where Lz. - z, - z, , and Az, . - = z, , , - z, . Note that If Az, , , = Az, , 
k k k-1 k+1 k+1 k k+1 k* 

the approximation collapses to the ordinary second-order centered 
approximation. At the surface, the vertical derivative Is diagnostic 
due to the surface boundary conditions; a one-sided difference Is used 

~ ^2 -Xi 

\3z/j^ Azj^ 

where Az^ Is the depth of the surface layer. At z^, again a one-sided 
scheme Is employed; ■ 



As mentioned earlier, the horizontal filter terms are evaluated 
using centered differences, l.e. 



The second derivatives In z needed for the vertical component of 
the sub-grid scale tendencies of u, v, 0^ and q are evaluated using the 
following centered difference; 
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ISzWsz/J^- ^^2 

where the has been evaluated at the half-grid levels, i.e. Az/2, 
3Az/2, •••. as calculated by Eq. C20) is not allowed to become 

greater than 0.85 of the stability criterion which for a centered-in- 
space forward-in-time scheme is K ^ V when the diffusion term is 

lagged in time. This limit is not reached in the runs described later. 

The basic time algorithm for all prognostic variables except r is 
the centered or leapfrog time step; thus. 



Because of numerical stability considerations, this scheme must be 
employed as follows : 

1. The advective and forcing terms remain unmodified and are evaluated 
at time level t. 

2. The sub-grid scale tendencies are lagged in time and are, there- 

(4) 

fore, evaluated at time level t-1. The stability limits on g 

( 2 ) 

and 3 discussed in Section 2.3.1 are modified by replacing At 


by 2At. Also the limiting of K discussed above is adjusted by 

X^ 

replacing At by 2 At. 


3. The convective scale tendencies are calculated similar to a for- 


ward method but are updated only every time steps and remain 
constant in between the updating (see Chapter 3). 

In order to control the time splitting associated with the leapfrog 
scheme a light time smoother is applied. The past value is adjusted 
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during updating as shown: 


x"-"- - x" + I [x^-^ + X^-^1 - 2X^] . 

where V Is 0.3 for and 0.1 for u, v, T^, q and c. 

2. 3. 2. 2 Differencing for rainwater 

The horizontal advectlve terms for r are approximated by upstream 

differences. If u , > 0, 

J 


u 



u 


l,j 



r 





while If u < 0 
1» J 


u 




At the boundaries, the r advectlve terms are diagnosed using the grad- 
ient just inside the boundary, for example at the west boundary, B; 

"b,i INb+i,j - • 

Note that depending on the direction of the wind, this scheme may be 
numerically unstable; however, this approximation does not feed back 
Into the model dynamics due to the lateral boundary conditions applied 
(see Section 2.3.3); it is calculated for diagnostic purposes only. 

The vertical advectlon terms for r also use the upstream scheme. 

The term in the rainwater equation which causes the most difficulty 
is the fall-velocity term, r V^)/9z. siLhce Is of order 10 m/s. 


3r 

9x 
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this term can become unstable for small hz. Therefore a seml-lmpllclt 
scheme Is used. Eq. (8) can be written as 


h. 

3t 


1_ 9<Pm^ V 

Pm 


f(r) , 


where f (r) Is the forcing term and Is equal to -V»Vr - w 5 — + conversion 

0 z 

rate + collection rate - evaporation rate of r. In seml-lmpllclt form 
using second-order upstream-ln-space and centered-ln-tlme differencing, 
we can write 


T+1 T T+1 T 

^k+1 ~ ^k+1 ~ 

2 At 





r 


T+1 

k+1 


^^m'^t^k *^k ^‘^m^t^k+l ^k+1 “ ^‘^m'^t^k ^k ^ 


or 


>T+1 


T ^ T 
^k ^ ’^k+l 


T+1 

^k+1 


At 


Az p 


m k+*4 


TT nT t+1 

^Pm\^k+1 ’^k+l 


^Pm'^t^k+l ’^k+l 


<Pm\>k ^ 


+ f(r)^ 2At 


/ 


, , at <PA>k 

^ Az |T 

Pm k+*5 


where p 


T 

m k-t^ 


Is the layer mean density at time step t; 


l.e. 



‘’mlfc + P-l 


m' k+1 


Note this Is an explicit equation which needs only an upper boundary 
condition to obtain a solution. The boundary condition Is ~ 0 
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for all time; therefore 



This scheme Is stable and non-damping. 

As with all seml-lmpllclt schemes, If the Courant-Frledrlchs-Lewy 
(C.F.L.) criterion Is violated the scheme remains stable but the phase 
speed Is slowed. This may be acceptable when the terms being handled 
seml-lmpllcltly are not of Interest; however. In this case the rainfall 
term Is of Interest and a reasonable phase speed Is necessary. To 
accomplish this, the above difference equation Is marched forward using 
10 small time steps for every mbdel time step; i.e. all the terms except 
the forcing terms are evaluated and stepped forward 10 times using a 
time step of At/10 while the forcing term is held constant. This small 
time step Is sufficient to keep the C.F.L. criterion from being violated. 

2.3.3 The Lateral Boundary Conditions 
Several authors (Davies, 1972; Sundstrtim, 1973; Ollger and 
SundstrOm, 1975; and others) have demonstrated that great care must be 
taken to Insure lateral boundary conditions which are mathematically 
well-posed for the hydrostatic primitive equations; however, to date no 
satisfactory solution to the problem has been demonstrated. Such well- 
posed conditions must not only be numerically stable but must be non- 
ref lectlve In order to be of practical value. 

Various practical formulations of boundary conditions have been 
used and discussed in the literature. Shapiro and O'Brien (1970) using 
the filtered equations and Williamson and Browning (1974) using the 
primitive equations designate a region as Inflow or outflow based on 
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whether the wind normal to the boundary is In or out of the grid domain. 
The boundary conditions are then based on this designation with large- 
scale quantities predominantly determining the conditions in inflow 
regions and model quantities predominantly determining the conditions 
In outflow regions. This scheme experiences difficulty when the large- 
scale winds specified at the boundary Indicate inflow, while the model- 
calculated winds one grid interval interior to the boundary indicate 
outflow. 

This dissertation develops a practical boundary technique which 
allows useful limited-area forecasts to be made in spite of the 
mathematically ill-posed nature of the problem. The set of time- 
dependent lateral boundary conditions presented allows large-scale waves 
to enter the limited-area forecast domain but does not allow exiting 
large- or small-scale waves to be reflected with sufficient amplitude 
to ruin the results. 

The boundary conditions consist of large-scale time-varying ten- 
dencies linearly combined with model-calculated tendencies. The large-, 
scale conditions may be either generated by a large-scale model for 
real-time forecasts or calculated from real-data analysis for research 
studies. In either case, it is usually necessary to interpolate the 
large-scale data in both time and space in order to obtain the boundary 
condition tendencies. 

The prediction of any dependent variable, x» can be written as 
follows: 


X^(I) - Xp(I) + W(I) 


3t 


At 


+ 1 - 


W(I) 


9X 

It 


Sis 


At 


I 


( 21 ) 
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where the subscripts n and p denote the "new" value after the boundary 
condition and the "previous" value at a former time. The s denotes the 
large-scale specified tendency of x* and m, the model-calculated ten- 
dency. The values of the weighting coefficients, W, are given below: 

{ 0.0 for I “ the boundary grid points , 

0.4 for I ■ the boundary - 1 grid point , 

0.7 for I - the boundary -'2 grid points , (22) 

0.9 for I = the boundary - 3 grid points , 

1.0 for I = all other interior grid points 

Thus, the value at the boundaries is completely specified by the large- 
scale Imposed value, while at a distance 4A from the boundaries, the 
variable value Is Identical to the model-calculated value. If 
is zero, Eqs. (21) and (22) are similar to the "sponge" boundary condi- 
tions used by Kessel and Wlnnlnghoff (1972). 

These equations are also similar to the boundary scheme used by 
the National Meteorological Center (NMC) limited-area fine-mesh model 
(LFM) (Cooley, 1973). The LFM uses a blend of the LFM and the six-layer 
hemispheric PE model tendencies on both u and v wind components, 
temperature and mass fields. The six-layer PE model tendencies are 
calculated from six-hourly PE history fields biquadratically inter- 
polated in space to the LFM grid points. The weighting is as follows: 
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W(I) 


f O.O for I ' 
0.33 for I 
0.67 for I 
1.0 for I ■ 


the boundary grid points » 

' the boundary - 1 grid point , 

' the boundary - 2 grid points , 
all other Interior grid points . 


The implications of Eqs. (21) and (22) can be Illustrated using the 
linear advection equation. 


3t 



(23) 


where C is a constant and x Is referred to as height throughout the 
remainder of this section. In the absence of external forcing, Eq. (21) 
may be written in continuous space as 

- W(>0 , (24) 

where W(x) is an analytic function which fits the distribution of W(I). 
Thus, Eq. (23) can now be written as 

|i + CK(*) ^ - 0 . (25) 

It can be seen that the effect of the boundary condition is to reduce 
the group velocity of a disturbance to zero as it approaches the boun- 
dary. Therefore, in effect, the boundary conditions as given by Eqs. 
(21) and (22) alters the basic advection equation such that the ill- 
posedness of the boundaries is no longer an issue. 



29 


Eq. (25) also implies that the disturbance’s wavelength must de- 
crease and the gradient of X must build up in the boundary region. This 
necessitates the use of a smoothing device in this region, A more com- 
plete discussion of the boundary conditions and their behavior is pre- 
sented in Perkey and Kreitzberg (1976). 
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3.0 THE CONVECTIVE PARAMETERIZATION SCHEME 

Previous studies, both observational (Rlehl and Malkus, 1961) and 
theoretical (Lilly, 1960; Charney and Ellassen, 1964; Ooyama, 1964, 1969 
Kuo, 1965; and others) have demonstrated that many tropical disturbances 
especially hurricanes, are driven by the effects of convection; latent 
heat release and vertical transport of heat, moisture and momentum. Due 
to the Importance of the convection it cannot be Ignored by numerical 
modelers; however, due to the large scale separation between the convec- 
tion and the tropical disturbance It Is not possible on today's com- 
puters to simulate both scales In detail. Thus, It Is necessary to 
parameterize the sub-grid scale effects of convection. 

The convective parameterization schemes to date can be divided Into 
three categories; the first two categorizations are based on the physi- 
cal mechanisms they emphasize and the third on the use of the scheme. 
Category 1 schemes emphasize the lateral mixing of cloud substance into 
the environment as the principal means of atmospheric sensible tempera- 
ture and vapor changes while Category 2 schemes emphasize the cumulus 
Induced environmental subsidence as the principal mechanism. Both types 
of schemes require a closure assumption to determine the cloud mass or 
percent areal coverage and Information about profiles of in-cloud 
variables such as temperature, moisture, mass, etc. Category 3 schemes 
or so-called convective adjustment schemes were developed mainly for 
general circulation or extended-range forecast models. These schemes 
emphasize the avoidance of unstable primitive equation integrations and 
the rapidity of the calculation. 
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Category 1 Includes schemes developed and/or used by Ooyama (1964, 
1969), Kuo (1965, 1974), Krlshnamurtl (1968), Rosenthal (1969, 1973), 
Krlshnamurtl and Moxlm (1971), Mathur (1974) among others. These 
schemes are all very similar to Kuo's (1965) scheme and thus, the basics 
of the category will be discussed by using this scheme as an example. 
Also, this scheme is chosen because it is probably the best known and 
the most used of the Category 1 schemes. The two necessary assumptions 
for Kuo's parameterization are; 1) the atmospheric sensible temperature 
and vapor changes are proportional to the cumulus cloud's excess temper- 
ature and vapor profiles and 2) the proportionality for mixing (the per- 

L I 

cent cloud cover) is proportional to the total moisture convergence in 

" i ' 

the column. As implemented by various researchers, the second assump- 
tion is often modified such as to use only the moisture convergence in 
the boundary layer. In the tropics, this yields a reasonable measure of 
the cumulus activity. 

To determine the cloud excess temperature and vapor profiles, Kuo 
originally assumed that the cloud temperature profile is that of the 
"moist adiabat through the condensation level of the representative 
surface layer" and that the vapor profile is the saturated specific 
humidity profile corresponding to that temperature. The cloud top is 
defined as the level where the cloud excess temperature becomes zero. 

It is also assumed that all the condensate falls out and is accumulated 
as surface precipitation. Kuo (1974) shows how the pseudo-adiabatic 
assumption can be modified to include entrainment effects on the excess 
temperature profile. This scheme, as first presented, made no attempt 
to account for the compressions! heating due to the cumulus induced 
environmental subsidence; however, Kuo (1974) argues that the subsidence 
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velocity is Included In the average large-scale vertical velocity and 
thus, the compresslonal heating is accounted for by this technique. 
Kuo bases the argument on his Eq. (2.9a) 
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and Eq. (2.9b) 
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where the c subscript refers to cloud or ascending air, the d subscript 
refers to cloud-free or descending air, is the convective latent heat 
release and u equals dp/dt. The fallacy in this argument is apparent in 
Eqs. (2.9) where 30/ 9p is used instead of 30^/ 3p and 30^/ 3p, respec- 
tively; i.e. Kuo has not taken proper account of the difference between 
30/ 3p in the cloud and in the subsiding environment and therefore, has 
not properly accounted for the influence of the subsidence heating on 
the vertical distribution of the heating. 

The closure assumption for determining the total cloud mass is 
based on the total large-scale moisture convergence in the column. That 
is to say the large-scale Integrated moisture convergence is balanced by 
the cloud moisture (vapor and liquid) requirements. This allows the 
calculation of the percent area of cloud cover which, with the indivi- 
dual cloud properties, yields the environmental sensible temperature and 
moisture changes Induced by the convection. The 1974 method differs 
significantly from the 1965 method in the treatment of the moisture 
flux (Kuo, 1965, Eqs. 3.8 and 3.11 vs. Kuo, 1974, Eq. 3.11); the 



Importance of which Is emphasized by Ceselskl (1973, discussion of 
Eq. 2.4). 
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Category 2 Is represented by the Arakawa scheme (Arakawa al. , 
1969; Ooyama, 1971; Arakawa and Schubert, 1974). As stated earlier, 
this scheme stresses the conservation of cumulus-scale mass; thus, the 
major mechanism for atmospheric sensible temperature and vapor change is 
through the cloud mass detralnment at the cloud top layer causing sub- 
sidence below that level. The closure assumption used by this scheme 
requires that the convection be In quasl-equilibrlum with the large- 
scale forcing. 

The Arakawa-Schubert cloud Is initiated by an updraft at the top 
of the mixed-layer; l.e. the convection parameterized by the scheme is 
boundary layer rooted convection. The cloud conserves static energy 
until the cloud excess virtual static energy vanishes; this level is 
defined as the cloud top. The effects of entrainment are assumed in- 
versely proportional to the radius, R, which is held constant with 
height. Therefore, clouds with larger R are less diluted by entrain- 
ment and tend to raise to greater heights than clouds with smaller R. 

The rainfall is parameterized simply as a fraction of the condensate 
formed by the cloud. The remaining condensate and all the cloud mass 
is then detrained in a narrow region at the cloud top level. This 
induces cooling and moistening at this level due to evaporation of cloud 
condensate and warming and drying at all other levels due to subsidence. 
The scheme as presented in Arakawa and Schubert (1974) includes a spec- 
trum of cloud radii. Because of entrainment, this also introduces a 
spectrum of detralnment levels. 


Ill 
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As stated earlier, the closure scheme assumes that the convection 

is in quasi-equilibrium with the large scale. This in essence means 

that the relaxation time of the atmosphere’s adjustment to convection 

is small compared to the time scale of the large-scale forcing. Arakawa 

and Schubert estimate that the time scale of the large-scale forcing is 
5 4 

"10 s while the adjustment is ~10 s. Thus, the adjustment process is 
in quasi-equilibrium with the large-scale and the cloud mass flux must 
be related to changes in the large-scale temperature and moisture 
fields, i.e. the changes in the atmosphere’s stability. 

Category 3 schemes include Manabe and Strickler (196.4) and Oliger 
et al. (1970) among others. These schemes were developed as a quick 
method to remove conditional instability when the atmosphere is moist. 
The scheme instantaneously adjusts the atmosphere such that the condi- 
tional instability is removed while conserving the vertical integral of 
the static or internal energy. 

Ceselski (1973) compared these three categories of parameterization 
and made the following observations; 

1. Vertical distribution of heating; 

(a) Category 1 (Kuo, 1965) heating occurs at all levels where the 
cloud temperature is warmer than the environment. The magni- 
tude of the heating is directly proportional to the cloud 
excess temperature. There is no cooling at any level. 

(b) Category 2 heating occurs at all levels except the detraln- 
ment level. The heating is proportional to the upward cloud 
mass flux and the environmental stability. At the detrain- 
ment level cooling occurs. Because of the spectrum of clouds, 
the temperature, change may be a combination of subsidence 
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warming from large clouds and evaporational cooling from 
small clouds. 

(c) Category 3 heating and cooling occur such as to stabilize the 
sounding. 

2. Vertical moisture distributions: 

(a) Category 1 (^uo, 1965) adds moisture to the environment at 
all levels between cloud base and cloud top. The moistening 
Is directly proportional to the cloud excess moisture. 

(b) Category 2 adds moisture at the detralnment level and drys at 
all levels below the detralnment level. As with the heating, 
the effect of the cloud spectrum is to cause some levels to 
experience moistening and drying due to different clouds, 

(c) Category 3 can add or subtract moisture at any level. 

As stated In the beginning of this section, there is no doubt about 
the importance of convection in the tropics where the convective latent 
heat release is the major source of energy for the large-scale distur- 
bances. In contrast, the major source of energy for mid-latitude dis- 
turbances is the large-scale baroclinicity . Thus, the Importance of 
the convection in the large-scale energetics is not as well defined. 
Tracton (1973) proposed that in some instances "cumulus convection plays 
a crucial role in the initiation of development through the release of 
latent heat in the vicinity of the cyclone center. In such cases, 
dynamical models that do not adequately simulate convective precipita- 
tion, especially as it might occur in an environment that is unsatu- 
rated, will fall to properly forecast the onset of development." He 
demonstrated this hypothesis by investigating 14 storms and noting that 
the convection near the low center apparently initiated development 
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earlier than would have occurred if only large-scale baroclinicity were 
operating. 

Independent of the convective effects on the mid-latitudes cyclone's 
energetics, it is necessary to include the convection before reasonable 
forecasts of quantitative precipitation can be made. As will be shown 
in Chapter 4, the convection strongly Influences the surface precipita- 
tion amounts both in the vicinity of the low and along the cold frontal 
band where squalls can occur. 

With the exception of the convective adjustment method, the other 
schemes were developed to parameterize surface-rooted convection. This 
restriction may be appropriate and non-limiting in tropical cyclones and 
hurricanes where most of the above authors have applied these schemes; 
however, this is a very limiting factor in extratropical cyclones. For 
example, these schemes are unable to account for the mid-level based 
Convection observed in the occluded and warm frontal regions of the 
cyclone (Kreltzberg and Brown, 1970). Because the frontal surface may 
be considered as a material boundary, this convection is independent of 
the surface variables. It is doubtful if even squall-line convection 
is surface rooted. Thus, a more general scheme is necessary for use in 
the mid-latitudes. 

Because latent heat released by convective clouds is an important 
source of energy in a conditionally unstable atmosphere in the mid- 
latitudes and because the vertical scale of the convection is the same 
order of magnitude as the vertical scale of cyclonic disturbances, it 
is necessary to parameterize properly the vertical distribution of heat- 
ing caused by convective latent heat release. The resulting meridional 
circulation of a large-scale system is different with low-level heating 
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than with upper-level heating. For example. If heating occurs In the 
upper troposphere, the Induced meridional circulation would consist of 
strong outflow in a shallow upper layer and weak inflow of air in a deep 
lower layer. However, if the heating occurs in the lower troposphere 
only, strong inflow would be expected in a shallow low layer with a weak 
outflow in a rather deep layer above. 

Although both of these types of circulation lead to conversion of 
potential energy into kinetic energy, their roles in cyclone development 
are quite different. The strong inflow in lower layers associated with 
low- and mid-level cloud depths enhances future convective activity and 
cyclone development. The deep convection associated with the mature 
stage of development provides the vertical energy transport to higher 
altitudes necessary for the maintenance of the cyclone kinetic energy. 
Thus, proper formulation of the cumulus cloud top and vertical distri- 
bution of sensible temperature and moisture changes is of significance 
in distinguishing different stages of weather system evolution. 

The above considerations lead to the conclusion that the para- 
meterization scheme should Include both horizontal diffusion of decay- 
ing cloud material and environmental changes due to cloud-induced sub- 
sidence. If either mechanism for cumulus-induced environmental change 
is neglected, the vertical distribution is altered. Improper calcula- 
tion of the subsidence as in Category 1 schemes yields convective 
induced changes which warm and moisten the environment throughout the 
cloud depth while neglect of the horizontal mixing of cloud debris as in 
Category 2 schemes overestimates the warming and drying at mid and lower 
cloud levels and the cooling and moistening at the cloud top level. 
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The vertical distribution of environmental moisture changes can be 
more important than the temperature changes in a conditionally unstable 
atmosphere because the state of the cloud base region humidity deter- 
mines whether convection will occur or not. In extreme cases, the dry 
return flow can actually dissipate the disturbance (Zipser, 1969). This 
necessitates a reasonable parameterization of the cloud efficiency in 
converting cloud water into rainwater. The re-evaporation of cloud 
water and its associated cooling is completely neglected by the pseudo- 
adiabatic assumption of Category 1 schemes and only superficially 
attended to by Category 2 schemes. The direct dynamical impact of con- 
vection depends on the net latent heating, including evaporational 
cooling as well as condensational heating. 

The above schemes are also designed for cyclonic-scale predictions 
where the convection is more nearly in equilibrium with the large-scale 
forcing. Therefore, both the Category 1 and 2 schemes use a closure 
which preserves this equilibrium. However, on the mesoscale the more 
transient nature of the convective systems with time scales ~3 h is of 
interest. Thus, a closure system that is more locally responsive is 
necessary. 


3,1 The Sequential Plume Convective Parameterization Scheme 

This scheme, as do the other schemes, requires a cloud model and a 
closure assumption. The cloud model is a one-dimensional Lagrangian 
cloud model. The closure scheme is based on a vertically integrated 
slice method. This is to say that the integrated cloud heating due to 
latent heat release is balanced by the environmental warming due to 
subsidence. This is a locally derived closure scheme which is based on 
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3 

the cumulus time scale (~10 s) being much shorter than the time scale 

4 

of the mesoscale (~10 s). 

As just stated, the cloud model Is a Lagranglan cumulus model with 
dynamics similar to those developed by Turner (1962), Squires and Turner 
(1962), Simpson et al. (1965), Weinstein and Davis (1968), Simpson and 
Wlggert (1969) and Danlelsen et al, (1972). The physical processes 
affecting the environment surrounding the cumulus cloud which have been 
Included In the model are environmental subsidence Induced by the build- 
ing cloud, sub-cloud evaporation of convective precipitation and hori- 
zontal mixing of the dlsslpdtlng cumulus cloud. 

The parameterization scheme uses finite-differences designed to 
conserve moisture and heat. The decrease In the vertical Integrated 
moisture (vapor and liquid) Is accounted for by the surface precipita- 
tion while the net Increase In temperature Is accounted for by the 
latent heating caused by the precipitation and net Increase In cloud 
water content. Small heat Imbalances arise In the model from releasing 
the heat of fusion and neglecting the heat required for melting. 

The most fundamental assumption needed to enable use of a 
Lagranglan cloud model for convective adjustment calculations is that 
the Initial buoyant parcel variable values are appropriate for each sub- 
sequent parcel; thus, by following the first parcel the vertical pro- 
files of the cloud variables become known, Warner (1970) and later 
Cotton (1975) maintain that one-dimensional models with inverse R 
entrainment cannot predict both cloud top and In-cloud quantities such 
as liquid water. Simpson (1971) refutes this argument stating that If a 
rainwater fallout mechanism Is Included as Is the case with the Environ- 

C'O . 

mental Meteorology Branch (EMB) one-dimensional model and the model used 
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in this par^eterlzatlon scheme, then both cloud top and In-cloud pro- 
perties can be predicted. It Is true that time-dependent cloud models 
with many precipitation categories do a better job of predicting the 
precipitation, but the computer time required prohibits their use In a 
parameterization scheme. Also, It should be noted that Lagranglan dyna- 
mics do not explicitly treat the flow of the environment around the 
rising parcle but rather parameterize the effects of the environment on 
the parcel through the use of virtual mass and/or aerodynamic drag 
coefficients. 

The convective calculations are performed sequentially in the 
following steps: 

1. Convective base is selected; starting at the level of maximum 
releasable instability, the resulting convection must exceed a 
minimum depth (usually 600 m). 

2. Initial updraft conditions are determined; vertical velocity 
impulse and radius are specified; temperature, pressure and 
moisture are equal to environmental values. 

3. Updraft thermodynamics, microphysics, vertical velocity and radius 
are calculated, including the effects of entrainment, water load- 
ing, virtual mass correction, etc. 

4. Percent cumulus coverage calculations are performed. 

5. Subsidence and its effects on the environment are evaluated. 

6. Precipitation, collection and sub-cloud evaporation are computed. 

7. Horizontal mixing of the mature cloud and subsided environmental 
profiles is performed using a height-dependent ratio of cloud 
mass to environmental mass. 
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8. The modified environmental sounding is re-examined for further 
convection. 

9. The final sounding is adjusted for hydrostatic balance and mass 
continuity. Calculations are made of the total convective contri- 
bution to changes in the environment. 


3.2 Convective Base Selection and Initial Updraft Conditions 

The sounding or vertical column in the primitive equation model is 
Interpolated to levels with vertical spacing, Az, of 500 m for the con- 
vective calculations. Having obtained the sounding variables on the 
convective vertical grid, the updraft base level is selected as the 
level with the maximum releasable instability. To calculate this quan- 
tity, a parcel at a given level is Imparted an initial vertical velocity, 
(usually 1 m/s). This parcel then rises dry or saturated adiabati- 
cally, depending upon whether it is unsaturated or saturated, without 
any effects of freezing, entrainment or water loading. Based on the new 
parcel virtual temperature, the parcel acceleration is calculated and 
a new vertical velocity is computed; 


dw 

2 


g Buoy dz 


(26) 


where 


Buoy 


T - T 
vp ve 

T 

ve 


(27) 


and the subscripts p and e refer to the parcel and environmental values 
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Notice that the left-hand side of Eq. (26) is the change in the 
parcel’s specific vertical kinetic energy. Thus, the right-hand side is 
also an energy form and can be shown to be related to the area between 
the environmental and parcel temperature curves on a skew-T, In p dia- 
gram (Haltiner and Martin, 1957, p. 63). Therefore, the layer releas- 
able instability, rij^, is 


KTOP 


ri. 


g Buoy^ Az , 


(28) 


for Buoy^ >0.0 and where KTOP is one grid point below the level where 
the vertical velocity drops to less than ^'^pmln specified as 

Wp^/2). This quantity is calculated for each layer. 

The updraft base level selection begins where rlj^ is the largest 

proceeding to levels with successively smaller rij^ until ri^^ becomes 

2 2 

less than a preset minimum (usually 12.5 m /s ) or until a significant 
cloud is built. A significant cloud has been arbitrarily defined as a 
cloud with depth greater than 600 m; this prevents the convective ad- 
justment from spinning its wheels building many shallow clouds which 
have little effect on the PE forecast. Normally the level of highest 
rij^ does support a significant cloud but at times additional factors 
such as entrainment kill the cloud with updraft base at this level while 
a parcel at a level with smaller rij^ leads to a significant cloud. 

After selection of the updraft base, the initial cloud parameters 
are determined. The updraft parcel quantities of virtual temperature, 
pressure and specific humidity are set equal to the corresponding 
environmental values. The parcel is given an initial vertical velocity 
impulse, w^^ (usually 1 m/s), and an initial updraft radius, (for the 
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experiments in Chapter 4, •= 3 km). This radius is necessary for the 

entrainment calculations. 

The arbitrary specification of initial cloud radius is the weak 
link in Lagrangian convection models and no satisfactory algorithm has 
been devised to parameterize this quantity. Tests have shown that the 
12 h modified sounding results are not sensitive to the initial vertical 
velocity and that the results are sensitive to the initial radius only 
when entrainment is a major physical factor limiting the convective 
growth. In the cooler extra tropica, entrainment and hence, the initial 
radius are not critically important, but in the warmer tropics where the 
convection is more dependent on entrainment, the specification of the 
initial radius is more crucial, especially to the timing of the convec- 
tion. 


3.3 Updraft Calculations 

The updraft quantities are calculated similarly to those in the 
Welnsteln-Davis Penn State Lagrangian cumulus model (Weinstein and 
Davis, 1968) accounting for thermal buoyancy, vertical momentum entrain- 
ment and mass entrainment with the effects of the buoyancy correction 
using a virtual mass coefficient and the water loading as suggested by 
Simpson and Wlggert (1969). The mlcrophyslcs are parameterized follow- 
ing Kessler (1969). 

The updraft equations are the following system of total differen- 
tial equations describing the cloud growth in Lagrangian particle 


dynamics form 
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where term 1 is the virtual mass parameter, term 2 the thermal buoyancy, 
term 3 the water loading, term 4 the momentum entrainment and term 5 the 
form drag. 
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(30) 


where term 1 is the dry adiabatic cooling; term 2 the moist adiabatic 

latent heating (L =0 when the parcel is unsaturated), term 3 the 
CS 

temperature entrainment, terra 4 the evaporational resaturation after 

entrainment, term 5 the heat of fusion during freezing (L^ = 0 when 

T > T^. or T < or q and q, = 0) and term 6 the sublimation 
vp fi vp ff cp hp 

(L = 0 except at the temperature where the vapor pressure is reduced 
s 

from saturation over water to saturation over ice), q denotes the 

s 

saturated specific humidity; and are defined below. 
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where is the initial freezing temperature (usually 258K) , the 

final freezing temperature (usually 248K) and t^ is the time constant 
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for the freezing rate. For these calculations, is assumed lO"^ m. 

Thus, the conversion from liquid to solid occurs linearly over the 
temperature range from to 
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Note that below cloud base q^ changes by entrainment only while above 

cloud base q equals q . 
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where term 1 is the condensation, term 2 the suspended water dilution 
by entrainment, term 3 the vapor dilution by entrainment, term 4 the 
conversion of suspended water to rainwater (only if q^^ ^ 
threshold value) and tern 5 the collection of suspended water by rain- 
water. , K„ and q are constants evaluated following Kessler (1969), 
K 2 is a constant of proportionality prescribed by Kessler which relates 
the diameter and terminal velocity of a drop. This constant has been 
modified for temperatures below freezing. The fall speeds of large 
frozen drops are very comparable to those of large liquid drops but the 
fall speeds of small frozen drops are approximately 30% of those for 
corresponding liquid drops. Because it can be asserted that for a 
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Marshall-Palmer distribution large values of rainwater have larger 


drops than small values of fall velocities can be parameterized 

in a very elementary fashion by setting K2 in the fall velocity cal- 
culation proportional to when < TF, l.e. if < TF, 

K_ = KlfO.2 + 0.12(q. -0.5)] where is as prescribed by Kessler. 

2 2 hp 2 
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(34) 


where term 1 is the precipitation which fell out of the parcel, term 2 
the precipitation dilution by entrainment, term 3 the conversion and 
term 4 the collection. Again has been modified from Kessler's 
suggested value to account for freezing in the fall velocity calcula- 
tions, if Tp < TF, K^ = K^[0.2 + 0. 12 (q^p-0. 5) ] where K^ is as suggested 
by Kessler. 
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(35) 


where q^^^ is the precipitation dropped or unloaded from the rising 
parcel. 


There is no compelling reason for using the constant R^ as the 

characteristic depth in term 1 of Eqs. (34) and (35) to calculate the 

precipitation dropped, q^^^* Other characteristic depths that suggest 

themselves are; the parcel depth, P,, which varies with w ; the up- 

d p 

draft radius R(z) which varies with w ® and the cube root of the 

P 

2 1/3 

volume, L , = (P,*R ) (Which is independent of w . Therefore, 

VOX d p 
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^hd ~ \7~ with or characteristic depth; 


*hd 


with Pj as characteristic depth; 


w 


*hd 


-j- with R(z) as characteristic depth. 


Because a rapid inverse response Of to will Increase the unload- 
ing as w' slows down and reduce the unloading as w increases, P, would 
p p d 

give greatest stability and R(z) the least stability. Also because the 
initial parcel depth can be defined only arbitrarily and because use of 
R^ has proven satisfactory, the -latter has been used. 
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where term 1 is the acceleration, term 2 the mass entrainment, term 3 
the density effect. The entrainment parameter, p, is 


, _ 2a 

^ R 


There is no unique requirement for varying the radius of the parcel 
during ascent in a Lagrangian or particle dynamics formulation. Never- 
theless, the updraft radius, R(z), is needed to calculate the entrain- 
ment rate, the area of subsidence and the mass in the updraft column. 
Therefore, the vertical distribution of subsidence and of environmental 
heating is significantly influenced by the modeling of R(z). Three 
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possibilities that could be used without violating any physics Include: 
(1) constant updraft radius (Arakawa and Schubert, 1974); (2) spherical 
updraft parcel whose radius would be determined from mass continuity and 
density calculations (Ooyama, 1971); (3) vertical mass flux continuity 
as in a steady-state updraft (Squires and Turner, 1962). 

The third alternative has been selected as being the n. ?t plausible 
and, hopefully, the most effective means of modeling R(z). t should be 
recognized that while this formulation of R(z) is the same ai>.^in a 
steady-state updraft, this model as a whole does not envision a steady- 
state flow or balance in any other sense. 

This set of equations is solved by stepping upward dry or moist 

adibatically, performing the microphysical computations, calculating 

the entrainment effects and then diagnosing the new vertical velocity 

and radius. Care must be exercised when the cloud excess temperature 

becomes negative while the parcel is still ascending rapidly. The 

deceleration may become large and cause R to expand too rapidly; thus, 

implying an excessive build-up of cloud mass. This anviling effect 

cannot be treated properly in a particle dynamics model so an upper 

2 

limit must be imposed; AR/Az is limited such that R cannot double in 
less than 250 m (except as controlled by the upper boundary condition). 
This artificial limit brings the velocity profile more nearly in line 
with that obtained by Holton (1973, Fig. 3) through explicit inclusion 
of a pressure perturbation term in a fully time-dependent one- 
dimensional model. 

To insure conservation of mass, an upper boundary condition is 
necessary. This condition is applied at a height which is defined 

as the height where w equals w . . For this discussion, let z, be at 

T> * T^Tn *1 n '1 
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the grid point below and be at the grid point above z^^p (z^Qp 

is determined by interpolation when is found by testing for 
Wp < . The value of Pjjjpr|.Qp is found by interoplation between 

levels 1 and 2; R.j,Qp is calculated using Eq. (36) with Az replaced by 
AZtop ■ ^XOP ~ ^1* condition is specified that the mass passing 

beyond z^;^, M^p^p, is 


”p|xOP “ ^•^mpTOP ^ *^TOP ^ 


where D ■ 250 m. The factor of two arises from considering the mass in 

an updraft which decelerates at a constant rate, d, to zero velocity in 

2 

a distance D while maintaining constant mass flux, p irR w = 

’ mp p 

2 

‘^mpT0P'^’^T0P''pT0P* these conditions. 
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where d = Wp^Qp/2D. Thus, 
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where 6t is the time taken for the parcel to travel from z.j,Qp to z^^p+D; 
therefore. 


6t - t(z^QP+D) - t(z^Qp); Wp^Qp 
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Therefore, 


d 6t; 5t 
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“pTOP “ ^^mpTOP ^OP ° ■ 


(39) 


Furthermore, w at level z_ is found from 
P 2 
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(40) 


unless the above yields a negative w _ in which case w „ E 0. 

pz pZ 

3.4 Percent Cloud Cover and Environmental Subsidence Calculation 


To this point in the calculations, only the initial updraft parcel 
computations have been made. In order to complete the cloud,, i^jis 
assumed that the vertical profiles described by the initial as^c^nding 
parcel determine the vertical profiles of the mature cloud, l.e. 
subsequent parcels which "fill-up” the cloud have the same history as 
the initial parcel. 

This assumption completes the cloud updraft calculations. The 
plume or cloud mass, M^, which passed through a given level is the cloud 
mass above that level less the mass entrained above that level; 


= 


/ ^TOP 

ffrPjjjp(z') R^(z’) - Mp(z’) u(z')j dz' . 


(41) 


Or in discrete form 


MpU) 


M (z+Az) + TTR (z+Az) p (z+Az) Az 

_E SE 

1 + y(z) Az 


( 42 ) 


The total cloud mass which passed upward through the base level, M.„, 

pB 

is given by Eq. (41) evaluated at z^. , 
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The total plume mass,iy<!^, In the grid-point volume is M _ times 

P pB 

the number of clouds, N , where 

P 


N 

P 



(43) 


and Apg is the total c],oud area at the base height and is the cloud 
base radius. Thus, 


^ — 5^ . (44) 

TI Rg 

Before the subsidence values can be calculated, a mass detralnment 
prof il’e must be determined. It is assumed that the detralnment profile 
of the decaying cloud is the same as the cloud mass profile of the 
mature cloud, M^Cz). This assumption envisions that the updraft stops 
when the first parcel reaches the top and that the cloud has the shape 
of R(z) which was described by the initial parcel. The convection 
builds and decays without any steady-state stage. This picture is more 
in line with the Thunderstorm Project (Byers and Braham, 1949) descrip- 
tion than with a steady giant thunderstorm or a hot tower with contin- 
uous upward flux. 

The subsidence calculations begins by defining the "fuel" mass, 

M^, which supplies the cloud; 


“f ' \ PmB 


(45) 


where is the total area, p^g is the moist density at the cloud base 
and 6zj is the depth of the fuel layer. Equating./^ and yields 
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a E ^ = 

PB 


mB 




M 


pB 


(46) 


where ot is the percent cloud cover at the base. 

pB 

Requiring that the computed subsidence yield 
scale vertical motion (the net vertical motion is 
Section 3.7) yields 


zero net cumulus 
accounted for in 


w =(l-a)w +aw =0 
c p cs p cu 


(47) 


or 


■ j .■ 


6z 

s 



9 


(48) 


where the subscripts s and u refer to subsidence and updraft values and 

the c subscript refers to cumulus values. The w’s are averaged over the 

spatial scale ailci lifetime of the cumulus cloud. If we define as 

s 

the total mass subsiding due to convection, we can write 


or 


.^g(z) 

e^#i<p(z) 


A^[l - ap(z)] PjjgCz) < 5 Zg(z) 

A,^ a (z) P„„(z) 6z (z) 
j. p mp u 


(49) 




P_„(z) 

ms 

P (z) 

mu 


( 50 ) 


Substituting from Eq. (44) and A p 6z and dividing both 

p S 8 IDS S 

A^ and A by A^ yields 
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6z 




(51) 


2 2 

where a (z) •• 1 - a (z) ■ 1 - a ^[R (z)/R ]. Note that once a _ is 
S p pB B pB 

determined then 6z (z) can be evaluated. 

s 

However, at this point neither a „ or 6z^ are known. To close the 

pii t 

problem, the assumption is made that the Integrated pressure change 
through the depth of the cloud due to the latent heat release is equal 
to the Integrated pressure change in the environment through the same 
depth due to the sensible temperature change Induced by subsidence. 

This assumption is consistent with the cloud envisioned earlier in 
that this is the same as requiring equality of the integrated subsided 
environmental and cloud densities; thus, the cloud mean buoyancy is 
zero. This is essentially an integrated slice method. With this 
closure assumption, the subsidence and the percent cloud cover are 
uniquely defined. 

Thus, we can write 



and 



P__(z) 

mp 


dz 


6p 


s 



P„„(z) dz 
ms 


(52) 


(53) 


Let p (z) be the initial environmental density profile before subsi- 
me 

dence, then p = p + Ap and p >• p + Ap • This leads to the 
’ ms me ms mp me ^mp 


constraint that 
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TOP 


Ap„^(z) dz 
mp 


f: 


TOP 


Ap (z) dz 
ms 


(54) 


Note that Ap^^ Is calculated during the parcel ascent; thus, in order to 
reduce the effect of large cumulus cloud water and precipitation values, 
the water loading effects on Ap^p are reduced such that the maximum 
liquid water used to calculate is 0.5 gm/kg. From the equation of 
state, neglecting the effects of the liquid water, we can write for any 
level z in the subsiding air 




AP, 


AT 


vs 


vs 


which reduces to 




AT, 


vs 


vs 


if Ap at a constant level z is assumed small during the subsidence pro- 
s 

cess. AT^^ can be approximated by (9T^/3z)^ '^^s* This neglects hori- 
zontal advectlon and the dlabatic temperature changes due to subsidence 
evaporation. Thus, as a first approximation 



APg(z) dz 



Pg(z)Ye(z) 


6Zg(z) 


dz 


or 


= - -2B 

IT 


£ 


'TOP 


Mp(z)Ye(z) 




dz 


(55) 
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where the moist densities have been replaced by the virtual density and 

(9T^/3z)^ has been replaced by Because R(z) and are of the 

same order of magnitude and a „ Is of order 0.1, as a first approxima- 

pB 

tlon 6p Is directly proportional to a This allows a quick itera- 
P pji 

tlve method for solving for and 6z^. Note 6z^ is the subsidence at 
the base level . 

A first guess of 6zl is made (usually 0.25 km), then a*„ and 6p' 

I pfi 8 

are calculated. 6z^ is then 


6z. 


<Sp. 


^ y ^ 

6p* • f 
s 


(56) 


This 6z, yields 6p nearly equal to 6p , thus indicating that the approx- 
X s p 

f ^TOP 

imatlon that J Ap (z) dz is proportional to a „ is justified. 

z„ ms pB 

Changes due to subsidence in variables representative of a layer 
must be calculated accurately if the balance of moisture and heat is to 
be acceptable. This requires more care than used in the calculations 
above. Note that Eq. (46) yields the subsidence of a parcel starting at 
level z; in actuality what is desired is the subsidence of a parcel 
arriving at z. This requires a Lagranglan calculation of an environ- 
mental variable, such as temperature, which will arrive at level z. 
Starting with Eq. (41) and evaluating at level + 6z^ (see Fig. 3), 
we can write 


M (z+6z ) 
P s 


Hp(z) 


I. 


z+6z 


[ttp (z')R‘ 
mp 


(z') - 


Mp(z')y(z’)] 


dz' , (57) 


The mass subsiding, M , at level z + 6z is 

s s 
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/ z+(SZs 

irlR^ - R^(z*)]p^^(z’) dz* , (58) 


where R^ is the radius 
Eq. (50) dictates that 


of the circle having an area equal to 
for each cloud 


A^. 


M^(z) 


mp 

P (z) 

ms 


M^(z) . 


Again calculating the subsidence, 6z , such that the net cumulus scale 

8 

vertical velocity is zero yields 


Mp(z) 


M (z+6z ) 
P s 


-z+6z 

/ ' 


[irp (z')R (z') - M (z*)y(z')] dz* 
mp br P 
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MpCz) 


^z+6Zgj- p (z’+6z ) 1 

^ p“^(z»4-5z;) 


/ z+6z 

' 


7T R (z’)p (z*) 

ms 


1 _ ^ 
p 

ms 


(z'+6Zg) 

(z'+6z ) 
s 


dz'. (59) 


Because 6z^ Is ~100 m, almost no error is made by neglecting the 

2 

second Integral in Eq. (59) and replacing the coefficient of ttR^ by 

p (z'). Thus, Eq. (59) becomes 
mp 


Mp(z) 


/•z+tSz 

f. 


4 


mp 


p(z')Mp(z') 


dz ' 


(59a) 


A quick method of solving the above integral equation to determine 6z^ 

has been developed and is discussed in Kreitzberg et al. (1974). 

As stated earlier, the new environmental values at level z are 

deteimilned by Interpolating to find the variable values at z + 6z (z) 

8 

and then lowering a parcel with these values to z. During this sub- 
sidence process the effect of evaporation of existing cloud water is 
evaluated. Thus, a new environmental sounding is constructed; this 
sounding will be mixed with the dissipating cloud sounding in propor- 
tion to their respective areal coverage. 

3.5 Calculation of the Dissipating Cloud Profiles and Surface 
Precipitation 

The cloud profile of temperature is assumed to be the same as that 
traced by the lead parcel. At this point in the calculation, the sus- 
pended water profile is that of an active cumulus, not that of a dissi- 
pating cumulus cloud. It is therefore imperative to consider the 
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changes to this profile which occurs between these stages. One of the 
major changes is due to the collection of the suspended water by pre- 
cipitation. 

In considering the process of cloud "aging," several questions must 
be answered: 

1. What fraction of the precipitation falls within the updraft and 
what portion falls outside the updraft? 

2. What fraction of the updraft has precipitation falling through it? 

3. What 1s the collection rate of the precipitation? 

4. What is the role of entrainment between the building and dissipat- 
ing stages? 

5. What effect does the continued autoconversion have on the final 
suspended water profile? 

6. What is the effect of other buoyant parcels on the moisture budget 
in the cloud? 

As a first approximation to the effects of continued entrainment, 
autoconversion and convection which could only be considered appro- 
priately in a time-dependent model, it is assumed that their net effect 
is zero, l.e. that the cooling and drying caused by the additional 
entrainment plus the removal of suspended water by autoconversion is 
canceled by the addition of suspended water by further convection. 

Both the fraction of precipitation falling outside the updraft 
core, and the fraction of updraft core without falling water, f^^* 

are considered to be functions of the distance below cloud top; 


f 

ro 


Cz) 


TOP 



0.2(z 
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(60a) 
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f 

CO 


(z) 


0.2(z 


TOP 



(60b) 


It should be noted that this parameterization does not consider one of 
the most Important pareimeters affecting these quantltltes, that is, wind 
shear. The collection of suspended water by the falling water is para- 
meterized after Kessler (1969). 

Evaporation of falling water into subsaturated updraft below the 
condensation level and into environment air helps to cool and moisten 
the sounding. This process occurs in three regions: 

1. In the area around the updraft with precipitation falling through 
it. 

2. In the unsaturated updraft core below the condensation level. 

3. In the subcloud region below the updraft. 

At this time, the effect of evaporation in region 1 is neglected. 

Before the total precipitation can be calculated the precipitation 

dropped, profile must be adjusted, was calculated following 

the actively rising lead parcel; thus, the rainwater dropped at each 

2 

level is the amount unloaded at that level by a parcel of mass ttR 

The cloud is "filled-up" by like parcels, each of which dropped an 

amount q, ,. Therefore, the total rainwater dropped at a given level is 
hd 

calculated as the amount dropped by one parcel at that level times the 
number of parcels which passed through the level; 


qv,.(z) = — 5 . 

T7R^(z)Azp^^(z) 

mp 

The rainfall concentration, q^. 


( 61 ) 


at any level is then calculated 


as 
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. j ^ r^TOP 




^^ I TJ^/-™*' 


or 


<Jr(z> 


q (z+Az) + i I 

mp z 


z+Az 




+ Aq „ (z' ) - Aq (z’ ) 
^coll ^evap 


R^(z')p„(z') dz* 

Tnp 


(62) 


At this point, a decision is made as to whether collection or evapora- 
tion should occur; if q^^ is greater than zero then collection occurs, 
otherwise evaporation occurs. Following Kessler, the collection equa- 
tion is 




Z^“^Z vrTT / l\ / 


dz , 


where the terminal velocity of the falling precipitation is given by 


V_ = K, q 
t 4 ^r 


.125 


and 


5.094, T > 273. 15K; 

( ' 

=> / 5.094 j 1.0-[0.8+0.12(q^-0.5)] 


273.15-T 


273.15-T 


flJ 


[, T ^273. 15K; 

t X p 


5.094[0.2+0.12(q -0.5)], T < T., . 

r p ti 
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Thus, the terminal velocity used for falling water Is a constant for 
temperatures below the Initial freezing temperature and above 273. 15K 
and has a linearly Interpolated value for temperatures In the Inter- 
mediate region. The above equations are appropriate If the precipitat- 
ing water Is assumed to fall through the layer In one large mass. 
However, In reality, the precipitating water from lower layers falls 
through the base before precipitating water from higher levels. Again 
as a first approximation to this process. It Is asstaned that the pre- 
cipitating water falls through the layer In n bundles of Aq^ (usually 
1 gm/kg) each. Therefore n Is defined as q^^/Aq^,. 

Integrating Eq. (62) and combining with the effects of 
and n yields 

(63) 

This equation determines the amount of water collected by the precipita- 
tion falling through a layer of thickness Az in bundles of Aq^. 

However, the question arises as to what happens In the top layer 
of the cloud where no precipitating water enters from above. If no 
collection is allowed to occur, then all of the suspended water in this 
layer must be evaporated and mixed into the environment, thus adding too 
much moisture and thereby over-cooling the environment. As an upper 
boundary condition, the collection In the layer from z,j,Qp “ Az to z.j.Qp 
Is 


1 - exp 


( 


„ * .875 . 

-K 2 K 2 Aq^ Az 


) 


(1-f 

CO 
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Still following Kessler, the basic parameterized rate of evapora- 
tion of falling water Into subsaturated air Is 


dt 


K,N 13/20 .13/20 


(64) 


where Is the niunber of raindrops per unit volume per unit diameter 
range and Is set equal to 10^ and Is the density of the environmental 

_3 

air and Is assumed constant at 1 kg m . It should be noted that the 
assumption that Is constant Introduces error Into the calculation 
because It is physically evident that the small drops are depleted by 
the evaporation, thereby reducing and also altering the Marshall- 
Palmer distribution. Also, the change in dq^/dt under conditions other 
than standard pressure and density have not been considered. This 
deletion becomes Important as the height of the cloud base Increases. 

Integration of the above equation yields 


Aq (z) - q (z) - 

^evap' 'jj.' ' 


qj(z 


.7/20 S^*^e~'^se^^‘^ 

/ ^ XT 


20/7 


(65) 


where 


3.5 X 10^ p _ 

o U e 


Evaporation also uses incremental precipitation fallout; thus, q^ 
in the above equation is replaced by Aq^ and the process is allowed to 
occur r) times or until the atmosphere below the cloud is saturated. 

This value is then added to the environmental specific humidity and the 
environmental virtual temperature is lowered to account for the evapo- 
ratlonal cooling. The atmosphere is then tested for supersaturation 
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and If necessary Incremental adjustment is accomplished. is the 

total incremental evaporation and is the amount subtracted from q^. 

The precipitation at cloud base, Pg, and the surface precipitation, 
Pg are calculated as follows: 


*^rB ^mpB 


and 


w 


**rS ^meS 


w 


where is equal to 1 gm cm 


-3 


3.6 Detralnment Calculations 


( 66 ) 


(67) 


The density weighted cloud profiles of virtual temperature, 
specific humidity and suspended water are horizontally (isobarically) 
mixed with the density weighted subsided environmental profiles. Be- 
cause the cloud radius varies with height the actual fraction of cloud 
area must be determined at each vertical level. Other than this slight 
complication, the mixing is straightforward. The suspended water is 
evaporated until the environment becomes saturated; any remaining water 
is considered to be environmental stratiform cloud. 

If Xp ®nd Xg any properties of the plume and subsided environ- 
ment, the mixed value, x^* is 

'Sn (1-a )p + a p 

P s p*"p 

o'i 


9 


( 68 ) 
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where 


a = a 


R^(z) 


pB 


3.7 Final Calculations 

The modified sounding Is then re-examined for convection by recal- 
culating the layer releasable Instability. If any rlj^ remains, the 
convective process Is repeated; If not, the total effects of the convec- 
tion are calculated and formed Into tendencies. It Is seldom that more 
than one cloud Is necessary to stabilize the sounding. 

Having deduced the temperature change produced by all the convec- 
tive processes except for the mean vertical motion associated with the 
heating, the pressure can be adjusted to return to hydrostatic balance. 
The mass can also be adjusted to reflect the net density change; this 
mass adjustment Is achieved by the net cumulus scale vertical motion, 

w , that has been Ignored In the convection calculations up to this 
s 

point (see Section 3.4). 

The value of w Is [see Kreltzberg et al. , 1974, Eq. (B-13)] 

g 


w 

s 


1 Ap 
Pm® 


(69) 


where Ap Is the pressure change at a level because of the hydrostatic 

readjustment. The value of t, the lifetime of the convective adjustment. 

Is taken to be 40 min. The sounding change due to vertical advectlon by 

w depends upon the vertical displacement, w t; so the effect of w Is 
s s s 
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Independent of tlie arbitrary selection of T. Notice that this expansion 
Is at constant pressure following the air motion . 

After the sounding has been modified by the lifting w^t, the 
departures of the final temperature, water vapor and cloud water pro- 
files from the Initial profiles give the cumulus changes to be Intro- 
duced Into the mesoscale primitive equation model. 

The upper boundary condition in the mesoscale model Includes 
specification of the local pressure tendency at the top of the model, 

(3p/3t) . The cumulus convection will contribute the effect of sub- 

s, 

grid scale motions to that upper boundary condition, 

's,z^ t 

This term accounts for the mass pushed out the top of the one-dimensional 
model in response to the cumulus heating. Of course, the result is a 
high pressure dome that will induce horizontal mass divergence in the 
three-dimensional model and thereby, cause the surface pressure to fall. 

At this point, the convective changes are smoothed using a 1-2-1 
horizontal smoother. This prevents the convection from exciting 2A 
horizontal waves in the PE model. Before these changes can be used in 
the PE model, they must be converted to tendencies by dividing by a 
characteristic time between calls to the convective routine. Thus, the 
question arises, how often does the convective routine need to be called 
in order to prevent the static instability from allowing destructive 
unstable amplification? The key quantity in determining the static 
Instability in a conditionally unstable atmosphere Is the degree of 
saturation, l.e., is the relative humidity 100% or not? It should be 
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noted that convection is actually initiated somewhat before the environ- 
ment is saturated. This is because of the initial vertical velocity 
impluse given to the cloud parcel. The question can be cast in the 
following form: 


< 1 9Rh\ 

RH at /pg “ RH at 


(71) 


when evaluated at the critical level of the sounding. 

Note that because the horizontal advection of RH is small compared 
to the vertical advection, we can expand Eq. (71) in the following form: 


1 aRH ^ w aRH 

RH at RH dt ■■ RH az 


(72) 


Utilizing the definition of relative humidity and saturation specific 
humidity, we have 




^ A ^ A J. ^ A 

— dq - de H dp 

q e s p 

^ s 


or, using the Clauslus-Clapeyron equation: 


(73) 
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The pressure term is neglected for this argument because it is small. 
Thus , 
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(74) 
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where dq/dt =• 0, dT /dt ■ -wF and 3T /3z = -v . Following the same pro- 

V V e 

cedure and noting that dT^/dt =* 0 for the convection changes, 


1 

3Rh\ 

i 92^ 

L 3^v\ 
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at / 

q at / 

r, -,.2 9t / 


conv 
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Thus, 
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The convective rates are calculated from 

6T 


3q 6q 

at “ V *At 
c 


9T 
\ 

at 


V *At 
c 


(76) 


where Is the convective routine calling frequency and At Is the 
PE time step. As an example, consider the following conditions: 


T 

V 

280K ; 

q = 

8,6 gm kg 

II 

5,9K km"^ ; 

-32 = 

3z 

-2,43 gm kg ^ km ^ 

w = 

13.8 cm s ^ ; 

r =» 

6K km ^ 

(ST = 

V 

0.45K ; 

6q = 

-0.55 gm kg~^ 

At = 

180 s 




Solving for convective frequency, we find - 14, Therefore, the con- 
vective routine must be called at least every 42 min In order to control 
the static Instability, may be more or less where, for example, w Is 
larger or smaller. 
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For the runs discussed in Chapter 4, the convective routine is 
called every 20 min. However, the changes are converted to tendencies 
prorating the changes over the subsequent 40 min and are added to the 
previous convective tendencies. Thus, at any one time, convective ten- 
dencies computed from conditions at two earlier times are being fed into 


the hydrostatic model. 
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4.0 A CASE STUDY 

During the period 1200 GMT 21 Febrjuary through 0600 (MT 
22 February 1971, an extratroplcal cyclone, initially centered near the 
southern tip of Texas, deepened from 996 mb to 993 mb and moved north- 
northeastward to the southwest corner of Missouri. As the cold front 
pushed eastward, the warm, moist, low-level air characterized by 
potential instability was lifted and severe convection broke out. The 
results were tornadoes, lightning, hail, surface winds as high as 50 to 
60 mph, and heavy rainfall amounts ahead of the front in east Texas, 
Louisiana, Mississippi and northward to Ohio. All told, the convective 
activity ahead of the storm resulted in 117 dead, nearly 1600 injured 
and over $20 million in damages (Decker, 1973). Because of its obvious 
severe convective activity, this storm was selected as an initial test 
case. 


4.1 Analysis and Initialization 

The temperature, wind field and surface pressure analyses were 
obtained from the National Center for Atmospheric Research (NCAR) Hemi- 
spheric General Circulation Real Data Model (GCM) . These analyses were 
in turn derived from the National Meteorological Center (NMC) tempera- 
ture and surface pressure analyses. The wind field was calculated geo- 
strophically from the hydrostatic pressure field and then adjusted to 
remove the vertically integrated mass divergence. This adjustment 
removes the Lamb wave as discussed by Washington and Baumhefner (1975). 
A complete description of the NCAR GCM is found in Oliger ^ al. (1970) 
and will not be discussed here except to point out that the version of 
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•the GCM used for these tests has a horizontal grid interval of 2.5° 

(-264 km) and a vertical grid interval of 3 km with a surface diagnostic 
layer of 1.5 km. 

To obtain the temperature and wind data at the fine-mesh grid 
points, the above fields were first horizontally interpolated using 
bicubic spline functions (Fulker, 1975) and then interpolated in the 
vertical. The temperature field was vertically interpolated linearly 
in z while the u and v components of the wind were interpolated using 
the exponential function shown in Eq. (77) between the surface and 
1.5 km and then linearly above. 

V(z) = ae^^ + c , (77) 

where V is either the u or v wind component and a, b and c are deter- 
mined such that V(0 km) = V(1.5 km) = V 2 and V(0.75 km) = 0.9V2* 

and V 2 are the appropriate GCM wind components at the surface and 
1.5 km. 

The surface pressure field is also interpolated to fine-mesh grid 
point data using bicubic spline functions. The pressure at all other 
levels is calculated hydrostaticlly from the interpolated grid point 
virtual temperature data. 

To obtain the initial specific humidity field, standard relative 
humidity rawinsonde data are analyzed using an optimal Interpolation 
technique following Gandin (1963). In general for any variable x» we 
can write 

N 

Xest ' Xg . 


(78) 
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where x *■ estimated grid point field of x» X the grid point 

cSt. g 

guess field, the station observed field, X^ the guess field bi- 

linearly interpolated to the station grid. N is the number of observa- 
tions within the range of influence of X; A is the weighting function 
which is chosen such as to minimize, in the least-square sense, the 
difference between the true value and the estimated value of X for a set 
of many observations, i.e. minimize “ ^est^^* Substituting 

X from Eq. (78) into this expression, differentiating with respect to 
and setting the results to zero yields the following matrix equation 


A 


where 


ilj = 


B 




^^^obs~^g^l^ ***^^^bs Xg^l^Xobs 

• • 

• * 2 


(79) 


B 


^^^obs ^g^l^^true ^g^l^ 
*^^^obs ^g^N^^true ^g^N^ 


Note if the x' field is the climatological field, i.e. x' = *^X^» then ip 
S S 

is the covariance matrix between the observations and climatological 

field. This matrix can be evaluated using several historical data sets 

or it can be modeled. The approach used here is to model the covariance 

-a d^ 


matrix as a decreasing function of distance, e , where a^ is chosen 

such as to reduce the weighting to 0.04 at d equal 500 km. The matrix 
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B cannot be evaluated using conventional data since the "true" values 

at the grid points are unknown. Therefore, it is assumed that it also 
-a^d^ 

follows the e decrease. 

Since in general the covariance matrix is invertable, we can then 
solve Eq. (79) for A, 

A = Bip~^ . 

Once A^ has been evaluated for each grid point then the estimated vari- 
able value is solved for using Eq. (78). 

In these experiments, the variable is the relative humidity and the 
guess field is not the climatological mean but is the Air Force Global 
Weather Central (AFGWC) humidity analysis. This perverts the Gandin 
statistical approach since it is not certain that the AFGWC analysis 
has the same statistics as the atmosphere. However, for the purposes of 
this research, this complication does not cause any difficulty. After 
the relative humidity has been determined by the above technique, the 
specific humidity is calculated using it and the interpolated tempera- 
ture field. 

The cloud and rain water fields are set to zero at all grid points. 
This is due to lack of data and does delay the onset of rain in the 
forecast. In addition to providing initial data, the GCM provides time- 
dependent boundary conditions for the fine-mesh model runs. Boundary 
conditions on u, v, T^, q and p^ are linearly interpolated in both time 
and space from the three-hour GCM history tapes. This procedure acts as 
a filter which allows only low frequency GCM information to be trans- 
mitted through the boundary. 
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4.2 Fine-Mesh Experiments 

The flne-mesh limited-area model has a 1.25® (~132 km) horizontal 
grid over the domain shown In Fig. 4. As stated above, 1200 QIT 
21 February 1971 was selected as the Initial data time. 

4.2.1 Experiment E-I 

The model for this experiment has a 1.5 km vertical grid Interval 

with a 25 m diagnostic boundary layer. Figure 5a shows the surface low 

located over the southern tip of Texas with a central pressure of 

996 mb. The Initial surface wind, temperature, specific humidity and 

relative humidity fields are shown In Fig. 5b through 5e. The cold 

front can be seen trailing south-southwest out of the low across 

-1 

Mexico; approximately between the 8.0 and 10.0 gm kg specific humidity 
Isollnes (Fig. 5d). As can be seen, there Is considerable low level 
advectlon of warm moist air over Louisiana, while at the same time, 
although not shown, there Is little moisture advectlon at 3 km. This 
differential moisture advectlon. In addition to an already potentially 
unstable atmosphere. Increases the possibility of severe convection 
occurring If a lifting source Is present. The lifting source Is, of 
course, the approaching cold front. 

Note that the winds are very weak (~1.5 m/s) and are nearly geo- 
stroplc In direction. In reality, the observed winds In the Gulf region 
were much stronger; for example, the observed wind at Lake Charles, 
Louisiana (30N, 93W) was ~8 m/s while the Initial conditions wind Is 
~1 m/s. Also, the Initialization has nearly eliminated the wind shift 
at the front; for example, the observed wind at Corpus Chrlstl, Texas 
(28N, 97W) was SSE at ~8 m/s, the Inltallzed wind Is SW at ~1 m/s. 
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Figure 4. Model domain. The dark solid line denotes the fine-mesh 
domain. 
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Figure 5. Flne-mesh initial conditions (1200 GMT). (a) Surface pres- 
sure Oiib), (b) Surface vector winds (m/s), Cc) Surface 
temperature C^C), (d) Surface specific humidity (gm/kg), 

(e) Surface relative humidity (%) . 
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while 2® west at Larado, Texas, the observed wind was NW at ~10 m/s and 
the initialized wind is NW at -1 m/s. 

The initialization also weakened the strength of the frontal 
temperature gradient; for example, the observed temperature at 
Brownsville,. Texas (southern most tip of Texas) was 23. 8C not ~19C as 
depicted in Fig. 5c; therefore, since the analyzed temperature in 
northern Texas is very nearly correct, the gradient is reduced by about 
5C in 10° latitude. 

During the subsequent 15 h, the low center moved first northeast- 
ward along the coast line and then northward to the southwest corner of 
Missouri (Fig. 6). The 0600 GMT surface pressure, wind, temperature, 
specific humidity and relative humidity fields are shown in Fig. 7. The 
observed positions and values of the central low for every 3 h between 
1200 GMT and 0600 GMT are shown in Fig. 8. The low during the first 6 h 
of the forecast tracked too far east before turning north; by 0600 GMT 
the position is very well predicted. It is speculated that this east- 
ward drift is partially due to the underestimation of the precipitation 
and therefore, the latent heat release north of the low during the first 
3 h to 6 h. This is suggested by the slightly farther eastward drift 
and less northerly movement of the low when the model was run dry, i.e. 
no latent heat release at all (Tracton et al, , 1975). 

The major cause of this underestimation of latent heat release is 
the weak vertical motions in the initial state after initialization. As 
discussed earlier, the GCM winds at all levels other than the surface 

i 

are initially geostrophic. These wind profiles are then adjusted to 
remove the vertically integrated mass divergence. pBecause this adjust- 
ment is small for geostrophic winds, the winds above the 1.5 km GCM 
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Fine-mesh forecast conditions (0600 GMT) . (a) Surface pres 

sure (mb), (b) Surface vector winds (m/s), (c) Surface 
temperature ("C), (d) Surface specific humidity (gm/kg), 

(e) Surface relative humidity (%). 
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boundary layer remain nearly geostrophlc and thus, contain very little 
convergence or divergence. The surface wind is geostrophlc in direction 
but is reduced in magnitude because of the GCM boundary layer drag 
formulation. Because of their weak nature, these winds contribute 
little to the convergence-divergence fields. The initial vertical velo- 
city fields at 4.5 km (~580 mb) have values between ±2 cm/s, while the 
equilibrium values attained by the model after 6 h of integration are on 
the order of ±10-15 cm/s. Thus, the model requires 3 h to 6 h to spin- 
up to equilibrium. 

Figure 9 shows the surface pressure at 30N, 92. 5W versus time and 
the central low pressure value versus time. The initial drop in pres- 
sure at 30N, 92. 5W is due to the approaching low center; however, from 
Fig. 9b it is apparent that the central pressure drops too rapidly 
initially. It also appears that there exists a small amplitude, low 
frequency gravity wave with a period ~6 h and amplitude ±1.5 mb. This 
wave is due to the initial imbalance between the mass and flow fields 
and is indicative of a problem with the initialization procedures. At 
this time, there is no satisfactory solution to this problem (Gerrity 
and McPherson, 1970). 

The wind, temperature and moisture patterns in Fig. 7 at 0600 GMT 
are indicative of the patterns during the 12 h to 15 h previous. The 
strong low level advectlon of warm moist air just ahead of the front has 
continued, thus creating and maintaining potential instability in the 
area. As stated earlier this, along with the frontal lifting in the 
region, creates large amounts of static instability which must be 
released by the convection. 



PRESSURE (mb) PRESSURE (mb) 



TIME (GMT) 
(a) 



TIME (GMT) 
(b) 


Figure 9. Surface pressure versus time. (a) Point 
30N, 92. 5W, (b) Central low pressure. 
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Hydrostatic instability was present initially over the Northern 
New Mexlco-Arlzona border and over Oklahoma, Arkansas and Louisiana; 
thereby initiating the convective parameterization routine. An example 
of this convection is shown in the next series of figures. The Initial 
sounding (1200 GMT) at 32. 5N, 92. 5W shows considerable potential and 
static instability (Fig. 10). Potential instability can be identified 
as a layer with decreasing static energy (CpT^ + gz + L^q); this condi- 
tion exists in the layer from 1,5 km to 4.5 km. Conditional instabi- 
lity can be identified as a region with decreasing saturated static 
energy (c T + gz + L q ) ; this condition also exists from 1.5 km to 

p V c^s 

4.5 km. Therefore, the potential instability can be realized in the 
fom of convection based at any level in this layer which becomes 
saturated or nearly saturated. The necessary condition used in the con- 
vective adjustment scheme is if a parcel which is given an initial 
vertical velocity impulse of w^^ can reach saturation before its verti- 
cal velocity decreases to then convection will occur based at this 

level. In this case a parcel beginning at 1.5 km with relative humidity 
of 95.4% reached saturation by 2.0 km (one convective grid interval) 

before w became less than w , * 
p pmin 

The resultant cloud is quite vigorous, with a depth of 9 km and a 
maximum vertical velocity of ~11.5 m/s (Fig. 11b). The increased 
temperature excess (Fig. 11a) at 8 km is due to the change of phase to 
ice. Note that in regions of acceleration the radius generally decreases 
while in regions of deceleration the radius expands (Fig. 11c). This 
would be the case exactly if there were no entrainment; however, since 
there is entrainment, the radius will expand even if there is no 
acceleration or deceleration. The cloud water or ice is a maximum in 
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Figure 10. Sounding for 32. 5N, 92. 5W at 1200 GMT. (a) Skew-T plot of 
temperature and dew-point, (b) Static energy (solid) and 
saturated static energy (dashed) profiles. 
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Figure 11. Cloud quantities for 32. 5N, 92. 5W at 1200 GMT. (a) Tempera- 
ture excess before and after subsidence, (b) Vertical velo- 
city, Cc) Cloud radius, (d) Liquid or solid water, 

(e) Induced environmental subsidence. 
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the lower cloud while the rain water or ice is a maximum in the upper 
cloud (Fig, lid). The sudden decrease in rain dropped at approximately 
8 km is due to the change from water to ice which has a somewhat, slower 
fall velocity than the liquid; thus, less falls out of the rising parcel. 
The cloud water after rain-out is re-evaporated into the environment or 
left as stratiform cloud water during mixing of the dissipated cloud 
debris. 

The percent cloud cover is 11.2%. Thus, to keep the net cloud 
vertical motion zero, the subsidence reached a magnitude of ~500 m in 
the upper regions of the cloud (Fig. lie). Figure 11a shows the temper- 
ature excess after subsidence. Note the integrated temperature excess 
is of order 0. This, as discussed in Chapter 3, yields a zero mean 
cloud buoyancy and an equal vertically integrated pressure change in the 
cloud and in the subsided environment. 

The rain from this cloud is 14.75 mm giving an area average precip- 
itation of 1.65 mm. This area average rain is assumed to fall over the 
next 20 min in the three-dimensional model. The cloud-induced tempera- 
ture and moisture changes cause the region at cloud base to warm and dry, 
leaving the potential instability but destroying the static instability 
by reducing the relative humidity at the cloud base. This result is 
depicted in both Figs. 12 and 13. It should be noted again that the 
sounding shown in Fig. 12 does not ever exist in the three-dimensional 
model, but rather the changes from the convection are Inserted over the 
next 40 min. Thus, if no other changes were occurring in the primitive 
equation model, the sounding in Fig. 12 would be the resultant sounding 
after 40 min of Integration; this, of course, is not the case as 
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(a) 


Figure 12. Cumulus modified sounding for 32. 5N, 92. 5W at 1200 GMT. 

Ca) Skew-T plot of temperature and dew-point, (b) Static 
energy (solid) and saturated static energy (dashed) 
profiles. 
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Figure 13. Convective induced changes. (a) Temperature change 
(b) Specific humidity change. 
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continued lifting is crewing new static instability as the old 
instability is being destroyed. 

After the initial out-break of convection, the convection withdrew 
to a region over Oklahoma and Northern Texas (Fig. 14). This band in- 
tensified and grew in southerly extent until 0000 GMT. During this time 
the band moved eastward to the Louisiana-Mississippi border. Between 
0000 GMT and 0600 GMT, the band continued to move eastward while de- 
creasing in intensity until by 0600 GMT the band had dissipated. 

Figure 15 shows the observed squall-line positions as determined 
from the NOAA Weather Service Hourly Radar Summaries. The positions 
shown are those of the major band, other lesser bands which formed and 
dissipated are not depicted. As can be seen, the model moved the squall 
eastward too rapidly. The predicted position at 0300 GMT corresponds 
approximately to the 0600 GMT observed position. 

The squall cloud tops were under predicted throughout the forecast 
period; however, they did exhibit the correct trends (see Table 1). The 

Table 1. Reported and Predicted Cloud Tops 


Time Reported Cloud Predicted Cloud 

GMT Tops Top 

Radar (Model) (1000*s ft) (1000* s ft) 


1145 (1200) 40-50 34.5 
1445 (1500) 32-40 20.6 
1745 (1800) 26-42 25.5 
2045 (2100) 43-54 27.2 
2345 (0000) 39-56 40.3 
0245 (0300) 31-50 27.2 
0545 (0600) 34-48 
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(b) 

A 

Figure 14. Convective precipitation rate [mm (10 s) ] at (a) 1800 GMT, 
(b) 2100 GMT, (c) 0000 GMT, (d) 0300 GMT, (e) 0600 GMT. 
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predicted convective precipitation rates and both the observed and pre- 
dicted cloud tops show the maximum intensity at about 0000 GMT. How- 
ever, the model dissipated the band far too early; observations show the 
band present at 0000 GMT, 23 February over Florida. The model’s under- 
prediction of cloud tops is due to the strong stable region above 
10.5 km in the initial GCM data. Recall the next data point in the GCM 
data was 13.5 km. Thus, there is a loss of resolution in determining 
the height of the tropopause. Actual data show the tropopause in the 
squall area located at ~200 mb not 250 mb as the initial data indicated. 

Stable precipitation did not begin until after 1700 GMT. This is 
contrary to the precipitation observations in that surface precipitation 
was being recorded in the Texas-Oklahoma area at 1200 GMT, i.e. at model 
start-up. This discrepancy can be explained as a fault of the initial 
data which did not contain any cloud or rain water. Therefore, the 
model required 3 h to form cloud water, convert it to rain water and 
then let the rain water fall to the ground. Figures 16, 17 and 18 show 
the accumulated convective and stable precipitation amounts at 1800, 

0000 and 0600 GMT. Figure 19 shows that although the precipitation rate 
is large over Louisiana and Mississippi, the model accumulation is 
generally less than 5 mm. This is partially because of the rapid east- 
ward movement of the squall and partially because of the lack of mois- 
ture over the Gulf and in the squall initiation region in Texas. 

Figure 19b shows the total accumulated precipiation at 0600 GMT 
while Fig. 20 depicts the observed 18 h surface precipitation amount. 

The forecast has correctly predicted the location of the maximum over 
Oklahoma but has underestimated its magnitude. The maximum over 
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Figure 17. Predicted Ca) convective and (b) stable accumulated pre- 
cipitation (mm) amounts at [>0000 GMT. 


I 

















130H 120H IlON lOOH 90H 80M 


(b) 

4 -1 

Figure 19. (a) Total precipitation rate [mm (10 s) ] at 0000 GMT. 

(b) Total 18 h accumulated' precipitation (mm) at 0600 GMT. 
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Observed 18 h period 21/12 to 22/06 GMT surface precipitation (mm). Contours are as labeled 
with central maximums as Indicated. 
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Kentucky and Tennessee Is completely missed by the forecast. Also 
the precipitation in the squall-line region has been underestimated. 

4.2.2 Experiment E-II 

As mentioned earlier, a narrow band of moisture from the squall- 
line formation region south along the Texas Gulf Coast is not present in 
the initial humidity field. Because of its narrow band structure, the 
standard rawinsonde network did not detect its presence; its existence 
is shown by a cloudy region in satellite pictures. To remedy this 
situation and test the sensitivity of the model to changes in relative 
humidity, bogus moisture soundings are inserted in the satellite cloudy 
area. The rawinsonde and bogus data are then analyzed as before. All 
other fields are identical to those in E-I. 

As can be seen from Figs.* 21 and 22, the effect of the bogus sound- 
ings is to moisten a band from Central Texas south to the Texas Coast. 
This Increased moisture caused an Increase in the precipltable water 
from 2 mm along the Texas-Oklahoma border to 26 mm just east of 
Brownsville, Texas (Fig. 23). 

This increased moisture initiated a larger stable precipitation 
4 -1 

rate [+2.75 mm(10 s) ] at 1500 GMT over the Texas-Oklahoma border. It 
also initiated convection off the Texas coast at the initial time step 
(Fig. 24a). The effect of this convection which was not present in E-I 
is to initiate a gravity wave with vertical velocity amplitude of 
±3 cm/s (Fig, 24b) and pressure change amplitude of ±0.2 mb (not shown). 
The wave could also be seen in the u and v wind components. The effect 
of this wave over the 3 h period was to warm and cool the mid-levels by 
±0.2C. Thus, the gravity waves excited by the convective adjustment 
scheme have small amplitudes and seem to cause little if any harm. 
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Figure 22. Relative humidity (%) difference fields, RH^ 
(a) Surface, (b) 4.5 km, (c) 7.5 km. 
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Figure 23. Initial precipitable water Onm) for (a) E-I, (b) E-II 
and (c) E-II minus E-I difference field. 
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(a) 
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Figure 24. 1500 GMT E-II minus E-I difference fields of (a) accumu- 

lated convective precipitation (mm) and (b) 4.5 km vertical 
velocity (cm/s). 
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By 0000 GMT, the accumulated stable and convective ivrecipitation 
amounts have Increased with the largest difference (9.7 mm) being the 
convective increase In the squall-line region (Fig. 25). The Increased 
precipitation and Its associated latent heat release has caused a lower- 
ing (-1 mb) of the surface pressure (Fig. 26a) and an increased 
(+6.5 cm/s) upward vertical velocity (Fig. 26b). The pressure change is 
due to the increased (1C) mid-level heating (Fig. 26d). The evaporation 
of rain water below cloud base has lowered (-2.6C) the surface tempera- 
ture (Fig. 26c). The advection of moist low-level air in the squall- 
line region has also been Increased; thus helping to self-propagate the 
squall line. This is shown in Fig. 27 by the Increase in the convective 
precipitation rate over northern Alabama, l.e. the squall line is still 
present in E-II at 0600 GMT. 

The increased humidity in the squall-line formation region and over 
the Gulf of Mexico has increased the squall-line precipitation by 10 mm 
while doing little to the stable and convective precipitation amounts in 
the region of the surface low (Fig. 28). The increased precipitation 
has deepened the low and increased the circulation which helped self- 
propagate the convection. 

4.2.3 Experiment E-III 

To test the effect of added vertical resolution in the low-level 
moist boundary layer, the model's vertical grid is changed to the 
following vertically expanding grid; 0.0, 0.025, 0.375, 0.75, 1.25, 2.0, 
3.0, 4.5, 6.0, 7.5, 9.0, 10.5, 12.0, 14.0, 16.5 km. This yields six 
grid points below 3 km Instead of the three used in E-I and E-II. The 
GCM initial data and time-dependent boundaries are interpolated to the 
above grid. The E-II enhanced relative humidity data is re-analyzed 


118 



(a) 
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Figure 25. 0000 GMT E-II minus E-I difference fields of (a) accumu- 

lated convective precipitation (mm) and (b) accumulated 
stable precipitation (mm) . 
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onto the above grid. Because the enhanced humidity data was used, 
comparisons are made between E-Il and E-III. 

The E-III squall-line formation region and the Louisiana Gulf Coast 
are slightly more moist (1 to 5 mm) while the Texas-Mexico Gulf Coast 
area is slightly drier (~3 mm) than the E-II case (Fig. 29). A cross- 
section at 35N (Fig. 30) shows that the additional boundary layer grid 
points have captured a low-level humidity maximum which was not detected 
in the E-II experiment. This added I6w-level humidity lowered by 1 km 
the convective bases over Oklahoma which resulted in Increased convec- 
tion early in the forecast (before 1800 GMT; see Fig. 31) which pumped 
more low-level moisture to mid- and upper levels; thereby helping to 
saturate the mid-levels and increase the vertical velocity. These two 
effects produced more stable precipitation in the enhanced convective 
region than occurred in E-II. 

The final E-III minus E-II precipitation difference fields 
(Fig. 32) exhibit increased precipitation amounts over the Oklahoma 
observed maximum; thus increasing the amount to within ~40% of the ob- 
served amount (Fig. 33). There was also increases (~2 to 4 mm) over the 
squall region. The major effect of the increased vertical resolution 
was to Increase the early convective moisture pumping and latent heat 
release in the region ahead of the approaching low. 
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(b) 

Figure 30. 1200 GMT relative humidity C%) cross-section at 35N lati- 

tude for (a) E-II and (b) E-III. The contour Interval Is 
10 and the shaded regions are areas of 100%. 
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Figure 32. 0600 GMT E-III minus E-II difference fields of (a) accumu- 

lated convective precipitation (mm), (b) accumulated stable 
precipitation (mm) and (c) total accumulated precipitation 














Figure 33. 0600 GMT E-III total accumulated precipitation (mm). 
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5,0 SUMMARY AND CONCLUSIONS 

5.1 Review of the Lateral Boundary Conditions 

Since present day computer economics do not allow unlimited calcu- 
lations, restrictions on the number of grid points and time steps are 
necessary when developing a numerical model. Numerical stability places 
still other restrictions on the relationship between length of time step 
and size of grid Interval. Therefore, if the phenomenon being studied 
requires a small grid interval and thus, short time steps, the model 
domain and/or the length of forecast must be limited. These restric- 
tions have made limited-area moddls a necessary tool for numerical 
investigation of intermediate- and small-scale phenomena. 

This limited-area restriction Introduces lateral boundaries to the 
model domain and thus, the problem of how to treat these boundaries. 

The ideal boundary would be one which seemed to not exist, i.e. one 
which did not allow the interior to know there were boundaries. Practi- 
cally, this means that the boundaries allow changes outside the limited- 
area domain to Influence the model interior while not letting the 
interior changes reflect at the boundaries. Mathematically, the re- 
quirements of the above conditions are not fully understood at this time. 

In spite of the mathematical uncertainties and problems, a practical 
solution to the choice of lateral boundaries is presented. This solution 
does not attack, but rather circumvents the mathematical requirements of 
a well-posed set of boundary conditions. The conditions mesh the large- 
scale Imposed tendencies and the model calculated tendencies. In addi- 
tion, a region of high wavenumber-selective damping is required. The 
major drawback to these conditions is that this region of increased 
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damping reduces the predictive usefulness of 15 to 20% of the grid and 
therefore, the usefulness of an equal amount of computer resources. 

These boundary conditions are investigated using analytic and 
numerical finite-difference techniques for advective waves and finite- 
difference methods for gravity waves. The results of these tests indi- 
cate that the boundary conditions break down long and intermediate length 
interior waves into short waves; thus, the need for the damping region to 
keep these waves from propagating back into the model interior. The 
conditions also allow energy on long and intermediate scales to enter 
the model domain with little degradation. 

Although not explicitly discussed, the results of Chapter 4 also 
demonstrate the utility of these conditions in that they allowed the 
forecast experiments to be conducted with little if any lateral boundary 
problems. The precipitation region over the Washington Coast (Figs. 19, 
28 and 33) is not a result of boundary problems as it might be inter- 
preted, but rather is a low pressure system beginning to enter the model 
domain through the western boundary as specified by the GCM boundary 
conditions. The circulation around the low is shown in the u and v 
components of the wind at 0600 GMT (Fig. 7a) and by the Increased mois- 
ture and warmer temperatures (Figs. 7b, c and d) off the Washington 
Coast. 


5.2 Review of the Convective Parameterization Scheme 

Again due to computer size and speed limitations, the scale of 
motions described by a numerical model must be limited. This necessi- 
tates the parameterization of scales that are too small to be resolved 
by the model grid. The purpose of any parameterization scheme is to 
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represent the effects of the phenomenon being parameterized on the 
phenomenon being calculated explicitly, i.e. the calculated phenomenon 
should react as if the parameterized phenomenon were also being calcu- 
lated explicitly. 

In this case, the scale of activity being parameterized is the con- 
vective scale. In mid-latitudes there is a wide variety of convection; 
surface rooted convection, mid-level based altostratus-embedded convec- 
tion, frontally-forced convection, etc. The parameterization scheme 
must be able to work in all situations, not for just one type of convec- 
tion. Previously developed parameterization schemes were constructed 
for tropical convection (low-level based and boundary layer forced) or 
for global general circulation models (long time-scale models). Neither 
of these types of schemes is applicable for short-range mid-latitude 
quantitative precipitation forecasts. 

The scheme presented is designed to have the versatility necessary 
in mid-latitudes and to be applicable for short-range forecasts. The 
results in Chapter 4 indicate that indeed the scheme is able to func- 
tion in the frontally-forced squall-line region, in the gently-rising 
altostratus region ahead of the approaching low center and, later in the 
forecast, in the over-riding region ahead of the warm front. The pre- 
dicted convective bases ranged from 1.7 km in the squall region to 
3.0 km ahead of the low. Also, the closure assumption is based on a 
shorter time scale which allows for shorter scale variations than pre- 
vious schemes. This is apparent in the one-shot convection off the 
Texas Coast in E-II. Instability was present initially with the 
enhanced humidity field and was therefore released, but since there was 
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no large-scale forcing the instability was not re-generated and there- 
fore neither was the convection. 

5.3 Review of the Case Study 

Three experiments are discussed, the first (E-I), or base case used 
a 1.5 km vertical grid Interval and humidity analysis based on standard 
rawinsonde observations. This experiment correctly predicted the posi- 
tion of the precipitation maximum over Oklahoma, but underpredicted its 
magnitude. A maximum over Kentucky-Tennessee was totally missed. The 
squall-line precipitation was also too light. The squall position was 
predicted to move eastward somewhat too rapidly. 

The predicted precipitation amount over Oklahoma was combined 
stable and convective in nature. This agrees with observations in that 
hourly precipitation records indicate steady rainfall rates over the 
forecast period with periods of heavy rainfall. In the squall region 
observations show rainfall periods lasting 1 h to 2 h with high rates. 
This type of precipitation was also indicated by the forecast. 

A narrow band of moisture extending from Central Texas south along 
the coast was suggested by satellite cloud observations. This band was 
too narrow to be observed by the conventional rawinsonde observations 
network. Case E-II attempted to enhance the initial moisture field to 
reflect this narrow band. The inclusion of the moisture band enhanced 
the squall-line precipitation while doing little to the Oklahoma maximum. 
The initiation and dissipation of the squall line as well as the squall- 
line precipitation amounts were affected by this narrow moisture band. 
Although the enhanced humidity field is not necessarily the "true" 
humidity field, it is a reasonable and possible field and thus, indicates 
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the sensitivity of short-range quantitative precipitation amounts to 
changes in moisture fields. 

The third experiment which used increased low-level vertical re- 
solution indicates that, even without more horizontal resolution and 
therefore without added observational costs, better short-range precipi- 
tation forecasts can be obtained. The Oklahoma maximum was increased by 
~25% so that it more nearly agreed with observations. Also some in- 
crease was noted in the convective region, which also improved the 
forecast . 

5.4 Suggestions for Future Research 

The most singular suggestion which peraeates this report is 
the need for better analysis and initialization techniques. The initial 
data from the GCM have been greatly maligned by this investigation; how- 
ever, in all fairness, the data have been somewhat misused and pushed 
beyond its intended purpose by the author. Several more sophisticated 
analysis and Initialization techniques are currently becoming available 
and are still being evaluated (for example; Schlatter, 1975; Shapiro and 
Hastings, 1973; Bleck, 1975). 

The convective adjustment needs to be further tested against real 
data such as the GATE and in the future, SESAME observations. These 
experiments have the increased observational time and space resolution 
necessary to yield data concerning convectlve-mesoscale interaction. 

Such questions as how do the parameterization scheme's thermodynamic 
energy, moisture and mass vertical transports agree with observed trans- 
ports need to be answered so improvements can be made. Also how this 
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scheme's transports compare with other scheme's transports needs to be 
examined . 

The original purpose for the development of this model was to con- 
struct a mesoscale model with a grid Interval on the order of 35 km over 
a domain ~1000 km square. Preliminary tests have been made on this grid 
but were beyond the scope of this report. Continued testing on 

this scale should be continued. 

Finally, the model's vertical coordinate should be changed to 
Include the effects of terrain. Work toward this goal Is also already 
underway. 
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APPENDIX 1 

It is reasonably transparent to understand the origin of the water 
loading in the pressure gradient terms of the momentum equations, 

Eqs. (1) and (2). It is true that the loading term is a smal^ correc- 
tion to the pressure; for example, writing the hydrostatic equation in 
finite form, we have 

Ap 

Al “ - Pm« • 

Using the definitions of and p as stated in Section 2.1 yields 

p^ = p(l + c + r) ; (Al.l) 

thus, p^ and p differ by an amount on the order of 0.1%. Therefore, 
even if the atmosphere were saturated and contained 1 gm/kg of liquid 
water through a depth of 10 km the "wet" surface pressure would differ 
from the "dry" surface pressure by only 1 mb. However, if we assume a 
grid interval of 100 km between the "wet" and "dry" sounding the pres- 
sure gradient is 1 mb/100 km as opposed to zero if we neglected the 
loading effects. Thus, the loading term is considered Important in 
evaluating the pressure gradient term. 

The loading terms in the temperature and pressure tendency equa- 
tions, Eqs. (3) and (4) are somewhat more obscure; thus, a derivation of 
these equations will follow. First, the pressure tendency equations; 
starting with the hydrostatic equation and differentiating with respect 
to t yields 
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(Al*2) 


where the loading effects are included so that the pressure at any 
level reflects the weight of the liquid and solid water above that 
level. Recalling Eq. (Al.l) and differentiating with respect to t 
yields 


. 9p(c+r) 

3t 9t 3t 


(A1.3) 


Neglecting the volume of the liquid water which Is approximately 
10~^ less than the volume of air, we can write the continuity equation 
as 


= _ ^.pv - ^ . 

3t 3z 


(A1.4) 


Combining Eqs. (A1.2), (A1.3) and (A1.4) and integrating with respect 
to z from the top, z^, down yields 


= 

3t 


- 8 





Adding and subtracting the quantity gj^ [^»p(c+r)^ + 3p(c+r)/9z] dz to 
the right-hand side yields 


-n 




m 


pd (c+r) 
dt 


■Jdz + I? - gP„wj^ 


lE 

3t 


(A1.5) 



138 


The inclusion of the spherical correction term In the continuity equa- 
tion yields Eq. (4). 

For the thermodynamic equation, we start with the first law for 
moist air 


pm 


dt p dt 


(A1.6) 


which neglects the heat capacity of the liquid water. Noting that 


^pm dt ^p^ ' dt ^p dt 


and substituting Eq. (A1.6) into the prognostic equation for T yields 


9T 


V 


at 



1 • 1 dp 

3z c ^ c p dt 
P 


or 


3T 

V 

3t 





1V+i(3£ 

P P\3t 


+ 



(A1.7) 


where again the hydrostatic equation with water loading has been used. 
Note this Is identical to Eq. (3) in Section 2.1 with the exception of 
the convective and small-scale effects. 

The loading terms in the diagnostic vertical velocity equation, 
Eq. (9), follow directly from the pressure tendency equation derived 
above . 


5 
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